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An Axisymmetric Analog Two-Layer Convective 
Heating Procedure with Application to the 
Evaluation of Space Shuttle Orbiter Wing 
Leading Edge and Windward Surface Heating 


1 Summary 

A numerical procedure for predicting the convective heating rate of hy- 
personic re-entry vehicles is described. The procedure, which is based on 
the axisymmetric analog, consists of obtaining the three-dimensional inviscid 
flowfield solution; then the surface streamlines and metrics are calculated 
using the inviscid velocity components on the surface; finally, an axisym- 
metric boundary layer code or approximate convective heating equations 
are used to evaluate heating rates. This approach yields heating predictions 
to general three-dimensional body shapes. 

The procedure has been applied to the prediction of the wing leading edge 
heating to the Space Shuttle Orbiter. The numerical results are compared 
with the results of heat transfer testing (OH66) of an 0.025 scale model of the 
Space Shuttle Orbiter configuration in the Calspan Hypersonic Shock Tunnel 
(HST) at Mach 10 and angles of attack of 30 and 40 degrees. Comparisons 
with STS-5 flight data at Mach 9.15 and angle of attack of 37.4 degrees and 
STS-2 flight data at Mach 12.86 and angle of attack of 39.7 degrees are also 
given. 


2 Introduction 

The calculation of aerodynamic heating on a three-dimensional body 
at hypersonic speeds is a challenging problem. Since wind tunnel testing 
cannot simulate the high temperature air environment of hypersonic flight, 
it is necessary to rely on computational fluid dynamic (CFD) flowfield code 
predictions. 

Numerical methods have been developed for solving thin-layer Navier- 
Stokes equations over complex three-dimensional geometries to calculate 
aerodynamic heating[l]. But even with the most advanced modern su- 
percomputer many hours of computer time are still required. A simpler 
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method to compute the viscous flow uses the w axisymmetric analog” for 
three-dimensional boundary layers developed by Cooke[2]. Following that 
approach, the general three-dimensional boundary layer equations are writ- 
ten in a streamline coordinate system, and the cross- flow is assumed to be 
small and can be neglected. This reduces the three-dimensional boundary 
layer equations to a form that is identical to those of axisymmetric flow, 
provided that the distance along a streamline is interpreted as the distance 
along an “equivalent body”, and the metric coefficient that describes the 
spreading of streamlines is interpreted as the radius of an equivalent body. 
This allows the existing axisymmetric boundary layer codes or approximate 
convective heating equations to be used to compute the approximate three- 
dimensional heating rate along a streamline. 

In order to apply the axisymmetric analog technique in computing ap- 
proximate heating of three-dimensional bodies, the inviscid surface stream- 
line paths and the metric coefficients need to be computed. DeJarnette and 
Davis[3] calculated the streamlines as the lines of steepest descent emanat- 
ing from the stagnation point. DeJarnette and Hamilton[4] developed a 
simple method for calculating streamlines from a known pressure distribu- 
tion. However, this approach has proven difficult to apply, unless the surface 
geometry and pressures can be described analytically. More success has been 
achieved in using the three-dimensional inviscid flowfield solution to com- 
pute surface streamlines and metric coefficients[5-9]. But, the disadvantage 
of this approach is that it requires more computer time than the engineering 
approximate methods described in references 3 and 4. The majority of the 
computer time is spent in obtaining the inviscid flowfield solution. 

Previous streamline codes[5-9] use three-dimensional flowfield predic- 
tions described either in a spherical or cylindrical coordinate system. How- 
ever, most recent CFD flowfield solutions use either Cartesian or general- 
ized coordinates. In this work, the Inviscid Equilibrium Computation in 
3-Dimensions (IEC3D), a three-dimensional, shock-capturing, inviscid CFD 
code[10], was used to compute the inviscid flowfield solution. A streamline 
code, which uses surface velocity components in Cartesian coordinates as 
input, has been developed to trace streamline paths and compute metric 
coefficients along the path. A boundary layer code, the Boundary Layer 
Integral Matrix Procedure(BLIMP)[ll], was used to evaluate the heating. 
Because of the failure of convergence while using the BLIMP code to eval- 
uate heating for the flight case, approximate convective heating equations 
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developed by Zoby et a/[ 12] were used to obtain heating rates. This code 
has been designated as the Axisymmetric Analog 2- Layer Convective Heat- 
ing (AA2LCH) code. A sample input is shown in Appendix A, and the 
listing of the program is provided in Appendix B. 

This document will provide a detailed description of the numerical pro- 
cedure used. The procedure is applied to an 0.025 scale Shuttle Orbiter at 
a wind tunnel condition of Mach 10 and angles of attack of 30 and 40 de- 
grees. The predicted heating rates are compared with experimental results, 
obtained in the OH66[21] test. The procedure is also applied to a trajectory 
point of the STS-5 flight at an altitude of 159,000 ft ,Mach number of 9.15 
and angle of attack of 37.4 degrees, and a trajectory point of the STS-2 
flight at an altitude of 179,920 ft, Mach number of 12.86 and angle of attack 
of 39.7 degrees. 


3 Inviscid Flowfield 

A three-dimensional, shock capturing, inviscid CFD code, IEC3D[10], 
was used to compute the inviscid flowfield solution. The IEC3D code is 
a general purpose three-dimensional Euler solution CFD code, which can 
compute inviscid flowfield solutions around general three-dimensional ge- 
ometries at a wide range of flight conditions and angles-of-attack. This 
code utilizes an upwind-biased, finite-volume, high order Total Variation 
Diminishing (TVD) scheme. Both Van Leer’s Monotone Upstream- centered 
Scheme for Conservation Laws (MUSCL) type flux- vector splitting[13-14] 
and Roe’s characteristic-based flux-difference splitting[15] are considered. 
An improved implicit solution algorithm called Lower-Upper Symmetric 
Gauss-Seidel (LU-SGS)[16] has also been incorporated. This code has been 
validated by wind tunnel and flight data[17]. 

4 Surface Streamlines 

Streamlines in a flow are defined as lines that, at any instant, are tangent 
to the velocity vector. The streamline equation in Cartesian coordinates can 
be written as 

dx 
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V 


dy 

dt 

dz 


( 2 ) 

(3) 


Starting at a initial location (xo» Vo, zo), equations (1) through (3) can be 
integrated to obtain streamline locations. 


5 Metric Coefficient 


A location (x,y,z) can be written in Cartesian coordinates as 

f = xi + yj + zk (4) 

Let (f , T), t) represent an orthogonal curvilinear coordinate system, where £ 
is in the direction of the streamline, r) is everywhere normal to the streamline 
on the body surface and r is normal to the body surface as shown in Figure 
1. On the body surface we have 

f=^ + j?e, (5) 

where and e v are unit vectors in the £ and T) direction, respectively. The 
unit tangent vector e v can be written as 


dr 



ffi+SJ+ 



dr 

Ihj 


( 6 ) 


where |JJ 
ordinates, 


is the metric coefficient. Since (£,f 7 ,r) represent orthogonal co- 
therefore we have 


^ r) — e T x (7) 

The velocity on the surface of the body can be written as 

V = 

= ui + vj + wk (8) 
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where u,u, w are Cartesian surface velocity components. Therefore 


u - v -r w r 

e t=v l +v J + v k 


( 9 ) 


The points on the body can be described as x = x(y y z). Therefore, the 
body surface can be described as 


F{x,y,z) = x - x(y,z) = 0 


( 10 ) 


The unit normal vector to the body surface can then be written as 


e T 


VF 

VF | 

r x i + r y j + r z k 


where 


and 



Equation (7) can now be written as 
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From Equations (6) and (14) we have 



dr 

dri 


1 f w u dx^ 

?J [ \v + v7h) 


(11) 

( 12 ) 

( 13 ) 


(14) 


(15) 
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and 


02 

dr 

35 


1 / v u dx\ 

{v + vd^J 


Let h = 


\WF 
I i 

= metric coefficient, then 
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h 

fw 
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and 


Combining Equations (17) and (18), we have 


h = 


VF 

( dz 

dy\ 

V 

Vdr, 

w d n ) 


(16) 


(17) 


(18) 


(19) 


Equation (19) can now be used to compute the metric coefficient, h, along 
the streamline. From DeJarnette[18], on the body surface 


d f dy\ _ d / dy\ _ d . . _ dv dy dv dz 

dt \dt]) dr] \dt) dr) U dydrj + dzdr) 

d fdz\ _ d ( dz\ _ d . . _ dwdy dwdz 

dt vdf?/ dt] \dt ) 8 t) W dy dr] + dz drj 

For a given h , initial values of |jj and || can be obtained from equations 
(17) and (18). Equations (20) and (21) can then be integrated to obtain 
5^ and || along the streamline. It can be shown that the heating rate is 
independent of the initial value of h. 

6 Computation of Surface Streamlines and Metrics 


( 20 ) 

( 21 ) 


A fourth order variable step Runge-Kutta integrator was used to inte- 
grate equations (1) through (3), (20) and (21). In order to integrate equa- 
tions (20) and (21) we need to know the values of and at each 

integration step. To compute JvfJ, we need the values of and also 
at each integration step. All of these values can be computed numerically 
from the coordinates and surface velocity components at each grid point. 
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Consider a master element (l and a mapping to an element fl e in the 
physical domain, as shown in Figure 2. Figure 3 shows the mapping from a 
master element to the entire physical domain. By using the finite element 
shape functions, the mapping from the (a, /?) domain to the physical (x, y, z) 
domain can be described as 

n 

x = ^x^a,/?) 

i = i 

n 

y = 

j= i 

n 

j=i 

Here (xj, yj,zj ) are the x, y, z coordinates of a local nodal point j in element 
fi e , n is the total number of nodal points of the element and ^ j is the shape 
function at nodal point j of the master element fi. Refering to Figure 2, the 
shape function at nodal points of each sub-element of the master element 
are Element I: 

#1 

*2 

*4 

Element II: 

*2 
*3 
*6 
*5 

Element III: 


*5 

1 

'o 

1 

t-H 

II 

(33) 

*6 

= a(l - P) 

(34) 

*9 

= a/3 

(35) 

*8 

= (1 - a)0 

(36) 


-(1 - a)0 

(29) 

— a/3 

(30) 

a(l +/?) 

(31) 

(1 ~ a )(l + P) 

(32) 


a/? 

(25) 

-(1 + a)/3 

(26) 

(l + a)(l + /J) 

(27) 

— (1 + f3)a 

(28) 


( 22 ) 

(23) 

(24) 
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Element IV: 


*4 

= -a(l-0) 

(37) 


= (l + a)(l-0) 

(38) 

*8 

= (1 + a)0 

(39) 

*7 

= — a0 

(40) 


The x coordinate of a point in the physical domain, which is the mapping 
of sub-element J, can be written as 


X = + X 2 $2 + X 5 tf 5 + 

= xiafi - x 2 (l + <*)0 + x 5 (l + a)(l + 0) - x 4 (l + 0)a (41) 

The derivatives and || are 
dx 

-fa = *10 -*20 + *5(1 + /?)- * 4 (l + /3) (42) 

dx 

— = x-iOt - x 2 (l + a) + x 5 (l + a) - x 4 a (43) 

At point 5 of element I (a = 0, 0 = 0) in the physical domain, the derivatives 
are 


dx 

~da s = Xs_X4 

(44) 

0X 

5a ls = - I2+Is 

(45) 


The x coordinate of a point in the physical domain, which is the mapping 
of sub-element II , can be written as 

X = * 2^2 + * 3*3 + * 6*6 + * 5*5 

= -x 2 (l - a)0 - x 3 a0 -(- x 6 a(l + 0 ) + *s(l - <*)(1 + 0) (46) 

The derivatives and |§ are 

*20 - *30 + *s(l + 0 ) - *s(l + 0 ) (47) 

-* 2 (1 ~ <*) - x 3 a + x 6 a + x 5 (l - a) (48) 


dx 

da 

dx 

00 
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At point 5 of element II (a = 0, (3 = 0) in the physical domain, the deriva- 
tives are 


= X 6 - x 5 (49) 

= -x 2 + x s (50) 

The x coordinate of a point in the physical domain which is the mapping of 
sub-element III, can be written as 

X = X 5 ^ 5 + X 6 $6 + 29^9 + £8^8 

= x 5 (l - a)(l - 13) + X 6 a(l - P) + x 9 a/3 + x 8 (l - a)/3 (51) 


dx . 
dp' 5 


The derivatives and || are 

= -x 5 (l-/?) + x 6 (l-/?) + x 9 /?-£ 8 /? (52) 

dx 

— = -x 5 (l - a) - x 6 a + x 9 a + x 8 (l - <*) (53) 

op 

At point 5 of element III (a = 0,/3 = 0) in the physical domain, the 
derivatives are 

|“| 5 = -Z5+X6 (54) 

= - 15+^8 (55) 

dp 

The x coordinate of a point in the physical domain which is the mapping of 
sub-element IV, can be written as 


X = X 4 '$! 4 + £5^5 + X8^8 + £7^7 


= -x 4 q(1 - P) + x 5 (l + q)(1 - P) 4- £8(1 + <*)P ~ x i a P 
The derivatives and |=| are 


dx 

da 

dx 

dp 


-x 4 (l - P) + x 5 (l - P) + x s p - 17 P 
x 4 q - £5(1 + a) 4- x 8 (l + a) - x 7 a 


(56) 
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(57) 

(58) 



At point 5 of element IV (a = 0,/3 = 0) in the physical domain, the deriva- 
tives are 


= - *4 + *5 (59) 

Sx 

^|S = -*5+ *8 (60) 


Therefore, the derivatives and || at point 5 in the physical domain are 




Similarly, we have 


-(^5 — I4 + Xq — X5 — £5 + *6 — x 4 + I5) 

^(l 6 - ^4) ( 61 ) 

\{~ x 2 + *5 ~ X 2 + *s - *5 + *8 - *5 + *g) 


^(*8 * 2 ) 


(62) 

dy _ 

da 

\(ye - Va) 

(63) 

dy 

dd 

2^8 - Dt) 

(64) 

dz 

da 

\( 2 g ~ 2 a) 

(65) 

dz 

da 

- *») 

(66) 


The partial derivatives of x,y,z with respect to a,/3 at all interior grid 
points can then be obtained from Eqs. (61) through (66). For points on the 
boundary, partial derivatives can be found as follows: Refering to Figure 
2, for points on the left boundary, only elements II and III exist. The 
derivatives and 0 at point 5 in the physical domain are 



dx 

da 


-(x 6 - *5 + *6 - *5) 

ar 6 - x s 

1 , 

-(-ar 2 + ar 5 - ar 5 + ar 8 ) 
^(*8-ar 2 ) 


(67) 


( 68 ) 
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For points on the right boundary, only elements I and IV exist. The deriva- 
tives J*j and || at point 5 in the physical domain are 


dx \ - 
da' 5 

i(l 5 - x 4 - x 4 + x 5 ) 


'Cal N 

Cn 

II II 

Xs - x 4 

i(-x 2 + x 5 - x 5 + x 8 ) 

(69) 

= 

^(*8- x 2 ) 

(70) 


For points on the upper boundary, only elements I and II exist. The deriva- 
tives and || at point 5 are 

dx. \ 

^|S = -(X5 ~ X 4 + X 6 ~ X 5 ) 

= \(x 6 -x 4 ) (71) 

dx 1 . 

qqU = -{-X 2 + X 5 -X 2 + X 5 ) 

= x 5 - x 2 (72) 


For points on the lower boundary, only elements III and IV exist. The 
derivatives and |j| at point 5 are 



1 , 



= 2 ^ — 15 + 16 — *4 + X 5 ) 



= i(x 6 -x 4 ) 

(73) 

dx 

1 , 


~d0^ 5 

= 2 ~ Xs + - x 5 + *s) 



= Xg — X 5 

(74) 


Similarly, we can calculate §£,§$,§- and fj| at all boundary points. By 
using similar finite element mappings, we can calculate |^,f 
and at all grid points. 

From the mapping, we have x — x(a,0),y = y(a,0) and z = z(a,0). 
Therefore 


dx 


dx 

da 


, dx 
da + dp dp 


(75) 


11 



dy 


(76) 


dz 


dy, dy 

^ia + ^dp 

dz dz 

d^ da+ Bp dl3 


(77) 


or 


dx 


• dx_ 
da 

dx - 

W 


dy 

dz 

= 

fe 

dz 

1 


U4 



33 J 



da 

d/3 


(78) 



Assuming there exists an inverse mapping, then we have a = a(x, y, z ) and 
f3 = /3(x, y, z ). Therefore 


, da , da , da , 
da - -~-dx + -x-dy + -r~dz 
dx dy dz 

d/3 = Tx dx + Ty dy + Tz dz 


(80) 

(81) 


or 



da 

dx 



da 

dy 

dp_ 

dy 


da 

dz 

dp 



dx 


dy 

■ 

dz 


(82) 


from Equation (79) 


da 

1 

d/3 

\T\ . 


(if) 2 + (If) 2 + (If ) 2 


^ dp J 

f dx dx . dy dy 


d0) 
dz dz > 


( dx dx j d\j_ dy _i_ dz_ dz_ \ 
ydadp ^ dadp da dp) 

(If ) 2 + (if) 2 + (If ) 2 . 


dx djt dz 

5cr da 
dx By dz 

m dp dp dp . 


dx 


dy 

dz 


(83) 
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where \ J\ is called the Jacobian of the transformation and is equal to the 
ratio of areas of element at physical plane to that at (a,/3) plane. | J\ must 
be greater than 0 to ensure the inverse mapping exists. 


\j\ = 


.dx. 


[® 2 + <|^) 2 + ( Hmg ) 2 + ( S ) 2 + (£) 2 ] 


,dz . 
d a' 


,dx , 

-w 


'fly, 

d p‘ 


. dz 

dp J 


dx dx dy dy dz dz 2 
dad P da dp dad /? 


(84) 


From equations (82) and (83) we have 


da _ J_ 
dx ~ |J| 
da _ J_ 

dy ~ \J\ 

da _ J_ 
dz ~ |J| 
dd _ J_ 
dx ~ \J\ 

dP = J_ 
dy ~ \J\ 
d£ _ J_ 
dz ~ \J\ 


' SO 2 

[— K— ) 2 

{ da [{ dp } 

&& 

1 

•g><g ) 2 


+ A 2 +(— ) 2 i 

+ ^ K dp’ J 

+ (lr) 2 + (l7) 2 ] 


da' 

+(ll) 2 + (fe) 2 ] 


<9a ' 
da ] 


+ (l!) 2 + (B 2 ] 


da 2 
. dz , 
da ' 


dx . dx dx 

Wfaxdp + 

dy . dx dx 
dp^dadp + 

dz . dx dx 

dp^Ip + 

dx dx dx 
fa { fadP + 
dy .dx dx 
d^ { d^cdp + 
dz dx dx 
da^^adP + 


dy dy 
da dp 
dy dy 
da dp + 
dy dy 
da dp 
dy dy 
da dp + 
dy dy 
da dp + 
dy dy 
da dp + 


gg« 85 > 

lg>i< 86 > 

gg>l< 87 > 

g|>l< 88 > 

gg>'< 89 > 

l| ),<90) 


Now, the derivatives J*, §f , f*, and C an be computed. 


dx 

dx da dx dp 

(91) 

dy 

da dy + dp dy 

dx 

dx da dx dp 

(92) 

dz 

da dz dp dz 

dv 

dv da dv dp 
da dy dp dy 

(93) 

dy 

dv 

dv da dv dp 

(94) 

dz 

da dz dp dz 

dw 

dw da dw dp 

(95) 

dy 

da dy dp dy 

dw 

dw da dw dp 

(96) 

dz 

da dz dp dz 
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The procedure for tracing streamlines and computing metrics is: 1) ob- 
tain and Ir at all grid points, 2) integrate equations (1) 

through (3), (20) and (21), and 3) use (19) to compute the metric. The 
streamline path does not always fall on the grid point, therefore interpo- 
lation is required to obtain velocity components and derivatives along the 
path in order to perform the numerical integration. A finite element based 
linear interpolation scheme is used. Consider a master element as shown in 
Figure 3. The shape function at nodal points of the master element are 


*1 

= (l-a)(l-/?) 

(97) 

*2 

II 

t— ‘ 
1 

3 

(98) 

*3 

= a/3 

(99) 

$4 

= (1 -a)/? 

( 100 ) 


Equations (22) through (24) describe the mapping from the master element 
to a physical cell. By knowing the x, y, z and Xj, y Zj values, equations (22) 
to (24) represent a system of non-linear equations of unknown a and /?, which 
can be solved by using either Newton’s method or successive approximation. 
The finite element shape function can then be used to interpolate surface 
velocity components and derivatives required for integration. 

7 Backward Streamline Tracing 

Surface streamlines emanate from the stagnation point and spread all 
over the body surface. It is very difficult to start a streamline path from 
the stagnation point and have it pass through a specific point on the body. 
Since we are interested in evaluating heating rates at specific locations on 
the leading edge, a procedure was developed to trace a streamline in a back- 
ward fashion, i.e. starting at a specific location and tracing the streamline 
toward the stagnation point. This can be done easily by reversing the sign 
of velocity components while integrating equations (1) through (3). Stream- 
line coordinates and the integration step size are saved at each integration 
step. The process can then be reversed at a point very close to the stag- 
nation point (the stagnation point can never be reached) and re-trace the 
streamline in a forward fashion and compute metric coefficients along the 
streamline path until the point of interest is reached. The streamline loca- 
tion, metric coefficient and pressure distribution along the streamline can 
then be used to evaluate the heating rate. 
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8 Heating Computation 


8.1 BLIMP 

The BLIMP codefll] determines the viscous flow field across the bound- 
ary layer around a blunt body for either an axisymmetric or planar shape. 
This code requires as input the distance from the stagnation point, met- 
ric coefficient and pressure ratio and computes the heating rate. Ideal gas, 
equilibrium flow and reacting gas chemistry may be used with the code. 
Convergence problems were encountered while using the code to evaluate 
wing leading edge heating rate for the flight case. The difficulty was associ- 
ated with the rapid pressure rise close to the wing leading edge. Therefore 
the BLIMP code was only used to obtain the heating rate for the wind 
tunnel cases. 

8.2 Approximate Convective-Heating Equations 

The approximate convective-heating equations developed by Zoby et 
at[ 12] were used to obtain heating rates for the flight cases. These heating 
rate relations, valid for both laminar and turbulent flow, have been shown to 
yield results which compare favorably with the more exact solution obtained 
from the BLIMP code for both nonreacting and reacting gas mixtures for 
either constant or variable entropy edge conditions. 

The laminar heating rate is computed from an equation which relates 
heating rate to the momentum thickness Reynolds number 

q L = 0.22(Ree, e )-\p m /p)(p m /p)p e u e (H aw - H w )(Pr w )~ 0 - 6 (101) 

where ()* quantities are computed using Eckert reference enthalpy relation, 
subscripts w , e and aw represent wall, edge and adiabatic wall respectively, 
and 9i , used to compute the momentum thickness Reynolds number, Re$, 
is given by the equation 

* 0l = 0M4( J* pyu e h 2 dS) 1/2 /(p et i e h) (102) 

where h is the metric coefficient, and s is the distance along a streamline. 
The Eckert’s reference enthalpy relation is given by 

H* = 0.5(H W + H e ) + 0.22(H aw - H e ) (103) 
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The boundary-layer thickness is given approximately by the equation 

(■ S/0) L = 5-55 (104) 

The turbulent heating is also computed from an equation that relates tur- 
bulent heating to momentum thickness Reynolds number 

q T = c 1 (Re e , e )- m (p m /p e )(fi*/n e ) m p e u e (H aw - S w )( Pt w )~ 0 - 4 (105) 


and 

8t = (c 2 f p'u e( p*rh«dS)V(P'*'h) (106) 

Jo 

m = 2/(N + 1) (107) 

ci = (1 /c 5 ) 2N/{n+1) [N/(N + 1)(N + 2)] m (108) 

c 2 = (1 + m)ci (109) 

c 3 = (1 + m) (110) 

c 4 = l/c 3 (111) 

c 5 = 2.2433 + 0.93iV (112) 

The value of N which is the exponent in the power law velocity profile 
relation, u/u e = (y/S) l ^ N , was computed from the expression 

N = 12.67 - 6.51og(i£e0, e ) + 1.21(log(/2e0, e )) 2 (113) 


which comes from a curve fit of axisymmetric nozzle-wall data[19], and 

(|)r = * + 1 + [( ^ 2 j^ + 1)(1 + 1.29(Pr w ) a333 ^)] (114) 

Since an inviscid solution is known, the boundary layer edge can be located 
through an iterative process of the momentum thickness equation, reference 
enthalpy equation and corresponding approximate ratio of boundary-layer 
thickness to momentum thickness. The inviscid properties at this location 
can then be interpolated to obtain the boundary-layer edge properties. 

9 Gas Properties 

The equilibrium air curvefits of Gupta et a/[ 20] were used to evaluate 
thermodynamic and transport properties of the equilibrium air. The prop- 
erties include enthalpy, total specific heat at constant pressure, compressibil- 
ity factor, viscosity, total thermal conductivity, and total Prandtl number. 
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These curvefits are valid from 500 degrees Kelvin to 30,000 degrees Kelvin 
over a pressure range of 10~ 4 to 10 2 atmospheres. 


10 Results and Discussion 

10.1 Wind Tunnel Case 

The AA2LCH code was used to calculate the leading edge heating to 
a 0.025 scale model Space Shuttle Orbiter configuration at a wind tunnel 
condition of Mach 10 and angles of attack of 30 and 40 degrees. The results 
were then compared with the OH66 heat transfer test data. 

Heat transfer test OH66 was conducted on a 0.025 scale semi-span Space 
Shuttle Orbiter model in September of 1976 in the Calspan 96-inch Hyper- 
sonic Shock Tunnel (HST) at a nominal Mach number of 10.2[21]. The 
right wing and a portion of the aft fuselage were deleted to allow a larger 
model size. The partial model was constructed of stainless steel. One of the 
test objectives was to obtain spanwise heat transfer rate distribution on the 
leading edge of the glove and wing. Figure 4 shows the layout of thin-film 
heat transfer gauges on the model. 

The first step in using the AA2LCH code is to obtain an inviscid flowfield 
solution. The quality of the CFD solution is highly influenced by the quality 
of the grid. In order to obtain good wing leading edge heating predictions, it 
is crucial to have good grid resolution in the shock-shock interaction region 
near the wing leading edge. Because the location of shock-shock interaction 
is not known a priori, the following procedure was adapted to ensure good 
grid resolution at the regions needed. An initial grid as shown in figure 5 was 
used to get approximate shock shape and location. The grid ha s dimensions 
of 81 X 61 x 75 in the /, J, K directions, respectively, where the I direction 
is from the nose toward tail, the J direction is from the surface of the body 
toward the outer boundary of the grid, and the K direction is wrapped 
around the body from top to bottom. An outer grid boundary adjustment 
program [22] was then used to push the outer boundary inward toward the 
shock. Figure 6 shows the resulting grid after the adjustment. Finally, the 
Multidimensional Self-Adaptive Grid Code (SAGE)[23] was used to adapt 
the grid based on the flowfield solution. Convergent solution is then obtained 
on the adaptive grid. 

After the inviscid flowfield solution was obtained, the streamline and 
metric coefficients from the stagnation region to each of the heat transfer 
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gauge location were computed by the backward streamline tracing proce- 
dure. Figure 7 shows the streamlines on the body surface. Heating rates 
were then computed by using the BLIMP code. Figure 8 shows the compar- 
ison between predicted and measured heating rates for the 30 degree angle 
of attack case. The result is plotted as q w versus x/Z, where x is the length 
along the x - axis . For the .025 scaled model, L is equal to 32.02 inches. 
Figure 9 shows the heating rate and pressure distribution along the leading 
edge of the model. In this double y plot, both q w and p/pstag are plotted 
against x/L, where p stag is the stagnation pressure. From this figure, one 
can see that the heating rate distribution follows closely the presure dis- 
tribution. Figure 10 shows a typical plot of pressure distribution along a 
leading edge streamline. In this plot pj p 3 tag along a streamline is plotted 
against the distance from the stagnation point. The rapid rise in pressure 
contributes to the higher heating rate in the leading edge area. From figure 
8, one can see that the predictions compare well with experimental data ex- 
cept at x/L = .6620 and x/L = .7703 locations where the predicted value is 
about 15% less than the measurements. Overall, the predictions match well 
with the measurements. Figure 11 shows the comparison between computed 
and measured heating rate for the 40 degree angle of attack case. Again, the 
predictions compare well with the experimental data except at x/L = .6530. 
A heating rate of more than 135 Bin/ ft 2 /s was measured while the predic- 
tion was about 62.5 Btu/ft 2 /$. The measured heating rate at this location 
is considered unreasonably high. 

10.2 STS-5 Flight Case 

Since the IEC3D code does not include chemical non-equilibrium effects, 
a point on the Space Shuttle STS-5 trajectory was chosen for comparison 
where the flow is expected to be at chemical equilibrium. The point is at 
an altitude of 159,000 feet. The Mach number is 9.15 and angle of attack 
is 37.4 degrees. Again, the IEC3D code and the grid adaption procedure 
described above were used to obtain the inviscid flowfield solution. The 
initial grid was the same as the wind tunnel case. During the Shuttle flight, 
the radiometer located near the leading edge panel 9 was used to measure 
heating rate in the wing shock-shock interaction area. The location viewed 
by the radiometer is at x = 1095 in, y = 256.488 in, and z = 291.94 in. The 
predicted heating rate using the approximate convective-heating equations 
of Ref. [12] is about 8.0 Btuj ft 2 /s at this location, and the heating rate 
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backed out from radiometer data is about 14.0Btu/ ft 2 /s[24]. Because of the 
large discrepancy between the prediction and measurement, a comparison 
of windward surface heating rate with flight data was made in order to 
determine the reason for this under-prediction. Thermocouple locations 
on the windward surface are shown in figure 12. The detailed description 
of thermocouple locations and heat transfer analysis of the STS-5 flight 
data can be found in Hartung and Throckmorton[25], The comparison of 
heating rates along the windward centerline is shown in figure 13. The results 
are plotted as q w versus x/L , where L is taken as 1280.8 inches. Good 
agreement with the flight data was obtained on the windward centerline. 
Lateral heating comparisons in planes normal to the centerline of the Shuttle 
Orbiter on the fuselage at x/L stations of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, and 
0.9 are shown in figures 14 through 21. The results are plotted as q w versus 
where y = 0 is on the centerline. In general, the prediction is within 15% 
of the measurement, except at three locations where the difference is larger. 
The comparison of heating rate on the surface of the wing is shown in figures 
22 through 27. The results are plotted as q w versus x/L at y equal to 184.8, 
233.6, 275.3, 322.7. 369.0 and 420.8 inches. There appears to be some scatter 
in the data. The trajectory point is at 1120 seconds after entry interface. 
By examining thermocouple data[25], one should realize that boundary-layer 
transition has occurred, and the flow has already turned turbulent on part 
of the wing. This contributes to the large scatter of the data. 

10,3 STS-2 Flight Case 

A point on the Space Shuttle STS-2 trajectory point was also chosen 
for comparison. The point is at an altitude of 179,920 feet, about 1090 
seconds after entry interface. The Mach number is 12.86 and the angle of 
attack is 39.7 degrees. At this altitude and flight condition, the equilibrium 
assumption is still valid. The grid used for this case has dimensions of 
151 x 61 x 110 in the I,J,K directions, respectively. The inviscid solution 
was obtained by applying the outer boundary adjustment procedure twice, 
the grid adaption procedure was not used in this case. 

The detailed description of the heat transfer analysis of the STS-2 flight 
data can be found in Reference [26]. Figure 28 shows the comparison of 
heating rate along the windward centerline. The results are plotted as q w 
versus x/L. Again, good agreement with the flight data was obtained on the 
windward centerline. Lateral heating comparison in planes normal to the 


19 



centerline of the Shuttle Orbiter on the fuselage at x /L stations of 0.1, 0.2, 
0.3, 0.4, 0.5, 0.6, 0.7, and 0.9 are shown in Figures 29 through 36 respectively. 
The results are plotted as q w versus y. Good agreement between prediction 
and measured flight data was obtained again. Figures 37 through 42 show 
the comparison of heating rate on the wing surface. The results are plotted 
as q w versus x/L at y equal to 184.8, 233.6, 275.3, 322.7, 369.0, and 420.8 
inches. Agreement between the prediction and flight data for the STS-2 case 
is better than that for the STS-5 case. 

The heating rate prediction at the location viewed by the leading edge 
panel 9 radiometer is 20.0 Btu/ ft 2 /s, and the heating rate backed out from 
radiometer data is 27.0 Btu/ ft 2 / s. This is certainly an improvement as 
compared with the STS-5 case, but there is still a 25% under- prediction. 

From the above comparisons, the under prediction of heating by the 
two-layer method in the wing leading edge shock-shock interaction area is 
probably caused by the deficiency of the inviscid flow solver to resolve the 
strong viscous effect in this region. Also, the flow in the vicinity of the wing 
leading edge is highly three-dimensional, and neglecting the cross flow in the 
boundary layer is probably another reason for the under prediction. 

11 Conclusions 

A procedure for evaluating surface heating rates for Orbiter-like vehicles 
at high angle of attack ha s been developed. The axisymmetric analog based 
procedure uses three-dimensional inviscid solutions in Cartesian coordinates 
to trace streamlines and compute metric coefficients, and uses either an 
axisymmetric boundary layer code or Zoby’s approximate convective-heating 
equations for heating rate computation. The procedure has been applied to 
the prediction of leading edge heating rate of an 0.025 scale shuttle orbiter at 
a wind tunnel condition of Mach 10 and angles of attack at 30 and 40 degrees. 
Good agreement was obtained between predictions and measurements. 

Heating rate predictions were also compared with flight-deduced data on 
an STS-5 trajectory point and an STS-2 trajectory point. Good agreement 
were obtained on the windward surface for both cases. However, the predic- 
tions on the wing leading edge shock-shock interaction area were lower than 
flight-deduced data. This discrepancy is probably caused by the inviscid 
solver being unable to resolve the strong viscous effect of shock-shock inter- 
action, and the neglect of cross flow in the boundary layer by the two-layer 
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method. 
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Figure 1: Streamline coordinates 
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Figure 2: Mapping from master element to a physical element 
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Figure 3: Mapping from master element to physical domain 
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Figure 5: Pitch plane initial grid 



Figure 6: Pitch plane grid after outer boundary adjustment 
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Figure 7: Surface streamlines 
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[eating Rate (Btu/ft**2/s) 


Orbiter Wing Leading Edge Heating Prediction 
OH66 Wind Tunnel Condition, Alpha = 30 deg. 



Figure 8: Comparison of predicted and measured heating rates, a = 30 deg. 


29 






Distance from Stagnation Point (ft) 


Figure 10: Typical pressure distribution along a leading edge streamline 
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Orbiter Wing Leading Edge Heating Prediction 
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Figure 11: Comparison of predicted a nd measured heating rates, a = 40 deg. 
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STS-5 Windward Centerline Heating Rate, Mach = 9.15, Alpha = 37.4 deg 
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Figure 13: Windward centerline heating distribution, STS-5 case, Mach = 
9.15, Alpha = 37.4 deg. 
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Figure 14: Circumferential heating distribution at x/L - .1, STS-5 case 
Mach = 9.15, Alpha = 37.4 deg. 
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Figure 15: Circumferential heating distribution at x/L = .2, STS-5 case, 
Mach = 9.15, Alpha = 37.4 deg. 
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Figure IS: Circumferential heating distribution at x/L = .5, STS-5 case. 
Mach = 9.15, Alpha = 37.4 deg. 



y(in) 

Figure 19: Circumferential heating distribution at x/L = .6, STS-5 case, 
Mach — 9.15, Alpha — 37.4 deg. 
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Figure 20: Circumferential heating distribution at x/ L = .7, STS-5 case, 
Mach = 9.15, Alpha = 37.4 deg. 
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Figure 21: Circumferential heating distribution at x/L = .9, STS-5 case, 
Mach = 9.15, Alpha = 37.4 deg. 
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Figure 23: Heating distribution on windward surface at y = 233.6 in., STS-5 
case, Mach = 9.15, Alpha = 37.4 deg. 
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Figure 24: Heating distribution 
case, Mac h = 9.15, Alpha = 37 
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Figure 26: Heating distribution on windward surface at y — 369.0 in., STS-5 
case, Mach - 9.15, Alpha = 37.4 deg. 
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Figure 27: Heating distribution on windward surface at y = 420.8 in., STS-5 
case, Mach = 9.15, Alpha = 37.4 deg. 
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STS-2 Windward Centerline Heating Distribution 



Figure 28: Heating distribution along windward centerline, STS-2 case, 
Mach = 12.86, Alpha = 39.7 deg. 
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Figure 31: Circumferential heating distribution at x/L = .3, STS-2 case 
Mach = 12.86, Alpha = 39.7 deg. 
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Figure 32: Circumferential heating distribution at x/L = .4, STS-2 case, 
Mach = 12.86, Alpha = 39.7 deg. 
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Figure 33: Circumferential heating distribution at x/L = .5, STS-2 case 
Mach = 12.86, Alpha = 39.7 (leg. 
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Figure 3-4: Circumferential heating distribution at x/L = .6, STS-2 case, 
Mach - 12.86, Alpha = 39.7 deg. 
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Figure 35: Circumferential heating distribution at x/L - .7, STS-2 case, 
M ach = 12.86, Alpha = 39.7 deg. 
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Figure 36: Circumferential heating distribution at x/L = .9, STS-2 case, 
Mach = 12.86, Alpha = 39.7 deg. 
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Figure 37: Heating distribution on windward surface at y = 184.8 in., STS-2 
case, Mach = 12. S6. Alpha = 39.7 deg. 
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Figure 38: Heating distribution on windward surface at y = 233.6 in., STS-2 
case, Mach = 12.86, Alpha = 39.7 deg. 


46 












0.96 0.98 1 1.02 1.04 1.06 1.08 

x/L 

Figure 41: Heating distribution on windward surface at y = 369.0 in., STS-2 
case, Mach — 12.86, Alpha = 39.7 deg . 
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Figure 42: Heating distribution on windward surface at y - 420.8 in., STS-2 
case, Mach = 12.86, Alpha = 39.7 deq. 





A Appendix A - Input Description 

Following is a brief description of input records for the AA2LCH code, Table 
A.l shows a sample input. 


record 

description 

i 

Grid file, must be written using 
FORTRAN UNFORMATTED I/O as follows: 
write(unit) idim, jdim, kdim 

write(unit) (((x(i j,k), i=l, idim) j=l jdim), k=l, kdim), 
(((y(i j,k),i=l, idim) j=l jdim), k=lkdim), 
(((z(ij,k),i=l,idim) j=l jdim),k=lkdim) 

2 

Density file, must be written using 
FORTRAN UNFORMATTED I/O as follows: 
write(unit) idim, jdim, kdim 

write(unit) (((rho(ij,k)j=l,idim) j=l jdim ),k=l, kdim) 

3 

X-component of velocity file, must be written using 
FORTRAN UNFORMATTED I/O as follows: 
write(unit) idim, jdim, kdim 

write(unit) (((u(i j,k), i=l, idim) j=l jdim),k=l, kdim) 

4 

Y-component of velocity file, must be written using 
FORTRAN UNFORMATTED I/O as follows: 
write(unit) idim, jdim, kdim 

write(unit) (((v(i j,k), i=l, idim) j=l jdim), k=lkdim) 

5 

Z-component of velocity file, must be written using 
FORTRAN UNFORMATTED I/O as follows: 
write(unit) idim, jdim, kdim 

write(unit) (((w(i j,k), i=l, idim) j=ljdim),k=l, kdim) 

6 

Pressure file, must be written using 
FORTRAN UNFORMATTED I/O as follows: 
write(uint) idim, jdim, kdim 

write(unit)(((p(ij,k),i=l,idim)j= 1 jdim ) ,k= 1 ,kdim ) 

7 

streamline output file, contains following: 
streamline coordinates, metric coefficient, 
distance from stagnation point and pressure 

8 

streamline outfile, suitable for BLIMP input, 
contains distance from stagnation point, metric 
coefficient and pressure to stagnation pressure ratio 
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record 

description 

9 

output surface grid file in PLOT3D format(unformatted) 

10 

an output PLOT3D q iUe(unformatted), contains 
surface heating rate, surface radiation equilibrium 
temperature, boundary layer thickness, momentum 
thickness, and skin friction 

11 

input file, contains wall temperature 

12 

a temporary file 

13 

an output file, contains location, distance from 
stagnation point, heating rate and skin friction 
along a streamline 

14 

output file, contains radiation equilibrium temperature 
along a streamline 

15 

output file, contains surface heating rate at grid point 

16 

outfile file, contains error message 

17 

grid scale 

18 

model scale 

19 

free stream pressure(psf) 

20 

free stream temperature (deg Rankin) 

21 

free stream velocity(fps) 

22 

free stream Mach number 

23 

equilibrium air flag 
= 1; equilibrium air 
= 0; perfect gas 

24 

turbulent heating flag 
= 1; turbulent heating computation 
= 0; laminar heating computation 

25 

radiation equilibrium temperature flag 
= 1; radiation equilibrium wall temperature 
= 0; constant wall temperature 

26 

surface emissivity 


wall temperature input flag 

= 1; wall temperature file is specified in record 11 
= 0; input initial wall temperature in record 28 

28 

initial wall temperature 

29 

angle of attack(deg.) 
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record \ description 










distance from stagnation point where streamline tracing 
starts(not used in backward streamline tracing) 


initial step size of integration 


maximum relative truncation error allowed in variable 
step Runge-Kutta integration 


minimum step size allowed in variable step Runge-Kutta 
integration 


number of integration steps to be taken(not used in 
backward tracing) 


number of streamlines to be generated(not used in 
backward tracing) 


initial streamline spreading angle(not used in 
backward tracing) 


backward streamline tracing flag 
= 1; yes 
= 0; no 


flag indicates whether initial location of streamline 
is on a grid point 

= 1; initial location is on a grid point 
= 0; initial location is input in the following record 


coordinates of initial location of streamline 


initial and final grid points for streamline tracing, 
always in the form of imin, imax, kmin, and kmax, where i 
is from nose to tail and k is wrapped around the body from 
top to bottom 


character indicates whether boundary layer parameter is 
available for interpolation along a streamline 
’y’ is available 
’n’ is not available 


a PLOT3D q file contains boundary thickness, momentum 
thickness, displacement thickness and Reynolds number 
per foot 


a PLOT3D q file contains REK values for k equals to 
.1, .25, .5, and 1. inch 
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' st s2_nunf . xyz 9 
9 sts2_nunf .ql' 

9 sts2_nunf ,q2' 

' sts2_nunf .q3' 

9 sts2_nunf .q4' 

9 sts2_nunf .p' 

9 sts2_ws.dat' 

9 sts2_wm.dat' 

9 sts2_basur .xyz' 

9 sts2_basur .q' 

'twall.dat' 

9 temp . dat 9 
'htrate.dat' 

9 rwall.dat' 

' sufheat .dat' 

'error.dat' 

'grid scale * 9 1. 

'model scale * 9 1 
'free stream pressure*' .84578 
'free stream temperature*' 470.88 
'free stream velocity*' 13674.5 
'free stream mach no. *' 12.85 
'equilibrium air flag =' 1 
'turbulent heating flag «' 0 
'radiation equilibrium flag' 1 
'surface emissivity' .907 
'wall temperature input flag' 0 
'wall temperature deg. R' 540. 

39.7 

10 . 

.01 

1.0e-4 

1.0e-5 

130 

1 

180. 

'backward integration flag*' 1 
'flag indicate grid point only' 1 
'starting location = ' 251.5 6.4 -60.6416 
' imin, imax,kmin,kmax' 140 140 110 110 
'transition parameter files' 'n' 

' /ceg/kwang/stream/roex3.q2' 

' /ceg/kwang/stream/roex3 .q4' 


Table 1: Sample AA2LCH input file 
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program aa21ch 
c 

c This program reads in an euler flowfield in Cartesian coordinate, 
c and trace the streamine and compute metric coefficient in a backward 
c or forward fashion 
c integrator is used 
c 

common /connty/ 
common /stagpt/ 
common /psidps/ 
common / noraxyz / 
common /grdxyz/ 

& 

& 

common /sfxiet/ 

& 

& 

common /dvdwyz/ 

& 

c 

c modification of 

c 

common /duxyz/ dudx(151, 110) , dudy (151, 110) ,dudz (151, 110) , 

& dvdx(151, 110) , dwdx (151, 110) 

c end of modification 

common /geom/ dxdy (151, 110) , dxdz (151, 110) 

common /surfpr/ vel, surf n, up, vp, wp, surf p, dxdzp, dxdyp, dxi, det 
common /sfbond/ msp (151) ,msr (151) , nsp (151) , nsr (151) 
c common /qxyz/ rho (151, 40, 110) ,u (151, 40, 110) ,v (151, 40, 110) , 

common /qxyz/ u(151,40,110) ,v(151,40,110) , 

& w ( 151, 40, 110) , p (151, 40, 110) 

common /save/ igsav, issav, isav, jsav, ksav, asav,bsav, gsav, 

& first , back, second, if atal 

common /freest/ pinf , tinf , rhoinf , vinf , vninf , cpinf , igas 
common /radeq/ itwall, wallt, irad, eps 
common /turbml/ iturb 
common /tempwl/ twall(1000) 
c 

c boundary and momentum thickness 
c 

common /bltrn/ deltal,thetl 
dimension y(50), us (151) 

dimension xb(1500) ,yb (1500) , zb (1500) , hb (1500) 
dimension q (151, 1, 110, 5) 

dimension ql 1(151, 1,110, 5) , ql2 (151,1,110,5) 
logical blank, first, back, second, fext 

character*80 gdf ile, rhof le, rhouf 1, rhovf 1, rhowf 1, pf ile, 

& sfile,mfile, twf ile, xf ile, qf ile, lqfl, lqf2, 

& ilf ile, out fie, rwf ile, shf ile, err fie 

character*30 dummy 
character*l flaura 
c 

c conversion factors 
c 

data cmin, fps, atm, den /2.54, 30.48, 2116.2, .01601846/ 
data rtok /1.8/ 
c 

c Stefan-Boltzman constant (w/cm**2/k) 
c 

data sigma /5.67e-12/ 
c 

c input surface grid locations 
c 

read (5,*) gdf ile 
read(5,*) rhofle 
read ( 5 , * ) rhouf 1 


:rom a given starting location. Third order Runge-Kutta 
,o perform numerical integration. 

il, jl,kl 

xstag,ystag, zstag 

psi (8) ,dpdxi (8) ,dpdeta (8) ,dpdzta (8) 

xn(151, 110) ,yn (151, 110) , zn (151, 110) ,xnp,ynp, znp 

xy z (151,40,110,3) , idim, jdim, kdim, 

iblank (151, 40, 110) , blank, isubs (2) , 

jsubs (2) ,ksubs (2) 

xixs (151,110) ,xiys (151, 110) ,xizs (151,110) , 
etxs (151,110) ,etys (151, 110) ,etzs (151, 110) , 
flgsm(151) ,flgsn(151) 

dvdy (151,110) ,dvdz (151, 110) ,dwdy (151, 110) , 
dwdz (151,110) 

8/16/93 



read (5,*) rhovfl 
read (5,*) rhowfl 
read(5,*) pfile 
read(5,*) sfile 
read{5,*) mfile 
read(5,*) xfile 
read(5,*) qfile 
read(5,*) twfile 
read (5, *) ilfile 
read(5,*) outfle 
read (5,*) rwfile 
read(5,*) shfile 
read(5,*) errfle 
read (5/ *) dummy, scaleg 
read (5,*) dummy, sea lem 
scale * scalem / scaleg 
c 

c input freest ream condition 
c 

read(5,*) dummy, pinf 
read (5,*) dummy, tinf 
read (5,*) dummy, vinf 
read (5,*) dummy, ami nf 
pinf = pinf / atm 
tinf = tinf / rtok 
vinf = vinf*fps 
read (5, *) dummy, igas 
c 

c input turbulent flag 

c 

read (5, *) dummy, iturb 
c 

c input radiation equilibrium flag 

c 

read(5,*) dummy, i rad 
c 

c input surface emissivity 

c 

read (5,*) dummy, eps 
c 

c input wall temperature flag 

c 

read (5, *) dummy, itwall 
if(itwall .eq. 0) then 
read (5,*) dummy, wa lit 
endif 
c 

c get free stream properties 
c 

iflag = 2 

if (igas .eq. 0) iflag = 0 

call enthalpy (pinf , hinf , tinf, iflag) 

if(ifatal .ne. 0) then 

write (6,*) ' Bad freestream condition encountered' 

stop 
endif 

write (6,*) ' hinf (curve fit) = r ,hinf 

call eqprop (pinf , tinf ,cpinf , rhoinf ,xmuinf ,xkinf ,prinf ,gaminf , 
& arminf , if lag) 

if(ifatal .ne. 0) then 

write (6,*) ' Bad freestream condition encountered' 
stop 
endif 
c 

c normalization factor for velocity (cm/sec) 
c 



vninf = sqrt (pinf *1 . 0133e6/rhoinf ) 


open (unit*7, f ile*gdf ile, form*' unformatted' , 

^status*' old' ) 

open (unit*8, f ile=rhof le, form*' unformatted' , 

^status*' old' ) 

open (unit=9, f ile=rhouf 1, form*' unformatted' , 
fistatus*' old' ) 

open (unit=10, f ile*rhovf 1, form*' unformatted' , 

Sstatus*' old' ) 

open (unit=ll, f ile=rhowf 1, form*' unformatted ' , 

&status=' old' ) 

open (unit*12, f ile*pfile/ form*' unformatted' , status*' old' ) 
open (unit=13, f ile=sf ile, form*' formatted ' , 

& status*' unknown' ) 

open <unit*14, f ile*mf ile, form*' formatted' , status*' unknown' , 
£err*9999) 

open (unit*18, file*twf ile, form*' formatted' , 

& status*' unknown' ) 

inquire (file=xf ile, exist =f ext ) 
if (f ext) then 

open (unit=16, file=xf ile, form*' unformatted' , status*' old' ) 
open (unit=17, f ile=qf ile, form*' unformatted' , status*' old' ) 
read ( 17 ) idim, jdim, kdim 
read (17) aminf , alp, re, time 

read(17) ( ( ( (q(i, j,k,nx) , i*l, idim) , j=l , jdim) , k=l,kdim) / 

& nx*l,5) 

else 

open (unit=16, f ile=xf ile, form*' unformatted' , status*' new' ) 
open (unit=17, file*qf ile, form*' unformatted' , status*' new' ) 
endif 

additional file 

open (unit*15, f ile=ilf ile, form*' unformatted' , status*' unknown' ) 
open (unit*24, f ile*outf le, form*' formatted' , status*' unknown' ) 
open (unit*25, f ile*rwf ile, form*' formatted' , status*' unknown' ) 
open (unit*30, f ile=shf ile, form*' formatted' , status*' unknown' ) 
open (unit=40, f ile=errf le, form*' formatted' , status*' unknown' ) 

idim = 0 
jdim = 0 
kdim = 0 
10 continue 

read (7) idim f jdim, kdim 

read (7 ) ( ( (xyz (i, j , k, 1) , i*l, idim) , j=l, jdim) , k=l, kdim) , 

& ( ( (xyz (i, j, k,2) , i=l, idim) , j=l, jdim) , k*l, kdim) , 

& ( ( (xyz (i,j,k,3),i*l, idim) , j=l, jdim) , k=l, kdim) 

multiply by scale factor 

do 20 j=l , jdim 
do 20 k=l, kdim 

do 20 i = l,idim 

xyz(i,j,k,l) = xyz(i,j,k,l) * scale 

xyz(i, j,k,2) = xyz(i,j,k,2) * scale 

xyz(i,j,k,3) * xyz(i,j,k,3) * scale 

continue 

modification of 1/11/94 

pi = atan2 (1 . 0, 0 . ) * 2. 
dangle = pi / kdim 
do 30 j * 1, jdim 
do 30 k = l,kdim 

angle = (k-l)*dangle 



xyz(l,j,k,2) = l.e-7 * sin(angle) * scale 
xyz(l,j,k,3) = l.e-7 * cos (angle) * scale 
30 continue 
c end of modification 
c 

c read rho 

c 

c read < 8 ) idim, j dim, kdim 

c read<8) { ( (rho (i, j, k) , i=l,idim) , j=l, jdim) ,k*l,kdim) 

c 

c read u 

c 

read { 9 ) idim, jdim, kdim 

read (9) ( ( (u (i, j,k) , i=l, idim) , j®l, jdim) , k=l,kdim) 

c 

c read v 

c 

read (10) idim, jdim, kdim 

read (10) (((v(i,j,k), i=l, idim) , j=l, jdim) , k=l,kdim) 

c 

c read w 

c 

read ( 1 1 ) idim, j dim, kdim 

read (11) (((w(i,j,k), i=l, idim) , j=l, jdim) , k=l,kdim) 

c 

c read p 

c 

read ( 12 ) idim, jdim, kdim 

read (12) (((p(i,j,k), i*l, idim) , j=l, jdim) , k=l,kdim) 

c 

c alp is angle of attack 
c 

read (5,*) alp 
pi = atan2 (0 . , -1 . ) 
dtr = pi / 180. 
alpr * alp * dtr 
cosa * cos (alpr) 
sina * sin (alpr) 
ksdim = kdim 

if (alp .It. 0.) ksdim = 1 
c 

c locate stagnation point (no yaw) 
c 

istag = 0 

do 160 i = l,idim 

us(i) *= u(i,l, ksdim) 

160 continue 
c 

c locate a pair grid points between which u changes sign, 
c 

do 170 i = idim, 2,-1 
iml = i - 1 

if (us (i) *us (iml) .It. 0.) then 
is = iml 
ispl * is + 1 
go to 180 

else if (usd) .eq. 0.) then 
istag = 1 
is = i 

ispl = is + 1 
go to 180 

else if (us (iml) .eq. 0.) then 
istag = 1 
is = iml 
ispl = is + 1 
go to 180 



endif 

170 continue 
c 

write (6,*) 9 Stagnation point cannot be located' 
stop 

180 continue 

write (6,*) ' is, ispl' , is, ispl 
c 

c find coordinates of stagnation point 
c 

xstag = 0* 
ystag * 0. 
zstag = 0 . 

if(istag .eq. 1) then 

xstag = xyz (ispl, 1, ksdim, 1) 
ystag = xyz (ispl, 1, ksdim, 2) 
zstag = xyz (ispl, 1, ksdim/ 3) 
else 
c 

c use linear interpolation to find coordinates of 
c stagnation point, 
c 

ds = sqrt ( (xyz (is, 1 , ksdim, 1) * xyz (ispl, 1, ksdim, 1) ) **2 + 

& (xyz (is, 1, ksdim, 3) - xyz (ispl, 1, ksdim, 3) ) **2) 

dsl * - us (is) * ds / (us (ispl) - us (is)) 
ds2 = ds - dsl 

write (6,*) ' ds, dsl, ds2' , ds, dsl, ds2 
c 

xstag = xyz (is, 1, ksdim, 1) + dsl/ds* 

& (xyz (ispl, 1, ksdim, 1) -xyz (is, 1, ksdim, 1) ) 

zstag = xyz (is, 1, ksdim, 3) + dsl/ds* 

& (xyz (ispl, 1, ksdim, 3) -xyz (is, 1, ksdim, 3) ) 

pstag = p (is, 1, ksdim) +dsl/ds* (p (ispl, 1, ksdim) -p (is, 1, ksdim) ) 
write (6,*) 9 Interpolated stagnation pressure : pstag 

endif 
c 

c compare with maximum pressure 
c 

do 190 i = l,idim 

if (p (i, 1,1) .gt. pstag) pstag = p(i,l,l) 

if (p (i, 1, kdim) .gt. pstag) pstag = p(i,l,kdim) 

190 continue 

write (6,*) ' Stagnation point is at', xstag,' 0 zstag 
write (6,*) ' Stagnation pressure is ', pstag 

c 

c find starting point location 
c 

read (5,*) dels 
if (dels .gt. dsl) then 
is s = is 
is = is - 1 
if (is .le. 0) then 
is = iss 
dels = dsl 
endif 

else if (dels .gt. ds2) then 
isspl = ispl 
ispl = ispl + 1 
if (ispl .gt. idim) then 
ispl = isspl 
dels = ds2 
endif 
endif 

if (alp .gt. 0.) then 

xstagl = xstag + dels/ds2 * (xyz (ispl, 1, ksdim, 1) -xstag) 
zstagl = zstag + dels/ds2 * (xyz (ispl, 1, ksdim, 3) -zstag) 



oo ooo ooo ooo ooo o ooo oooooo 


ystagl * 0 . 

xstag2 = xstag + dels/dsl * (xyz (is, 1, ksdim, 1) -xstag) 

zstag2 = zstag + dels/dsl * (xyz (is, 1, ksdim, 3) -zstag) 

ystag2 = 0 . 
else 

xstagl * xstag + dels/dsl * (xyz (is , 1, ksdim, 1) -xstag) 

zstagl * zstag + dels/dsl * (xyz (is, 1, ksdim, 3) -zstag) 

ystagl * 0 . 

xstag2 * xstag + dels/ds2 * (xyz (ispl, 1, ksdim, 1) -xstag) 

zstag2 = zstag + dels/ds2 * (xyz (ispl/ 1, ksdim, 3) -zstag) 

ystag2 = 0 . 
endif 

if (is .ne. iss) is * iss 
if(isspl .ne. ispl) ispl * isspl 
c 

c set up boundary values for each element 
c 

ms per = 0 
nsper = 0 
msl * 1 
nsl * 1 
mdim = idim 
ndim = kdim 

call bounds (mdim, ndim, msl, mdim, nsl, ndim, msper, nsper) 
c 

c compute surface metrics 
c 

call surmet (mdim, ndim, ms 1 , mdim, nsl, ndim, ms per, nsper) 
c 

n is number of functions to be integrated 
n - 5 

modifiation of 8/16/93 
n * 6 

h is initial step size 

read (5,*) saveh 
hsave = saveh 

back = .false. 

emax is maximum relative truncation error 
read(5,*) emax 

hmin is minimum step size allowed 
read (5,*) hmin 

nstep is the number of steps of integration 
read (5,*) nstep 

nsl is total number of streamlines to be generated 
read (5,*) nsl 

angi is initial spreading angle between streamlines 
c 

read (5,*) angi 
c 

c read in backward streamline tracing flag 
c 



read (5, *) dummy , iback 


c 

c 

c 

c 


c 


c 


c 


c 


itype = 1 : multiple streamline, imin, imax, kmin,kmax required 
ityp>e = 0 : single streamline 


read (5, *) 
read (5, *) 
read (5, *) 
read(5, *) 
if (f laura 
read (5, 
read (5, 


dummy, itype 
dummy, xst, yst, zst 
dummy, ibmin, ibmax, kbmin, kbmax 
dummy, f laura 

.eq.'y' .or. f laura .eq. ' Y' ) then 
*) Iqfl 
*) lqf 2 


open (unit=21, f ile=lqf 1, form=' unformatted' , status=' old' ) 
open (unit*22, f ile«lqf2, form=' unformatted' , status^' old' ) 


read (21 ) ildim, jldim, kldim 
read (21) duml , dum2 , dum3 , dum4 

read(21) ( ( (qll (i, l,k,nx) , i=l, ildim) ,k=l, kldim) ,nx=l, 5) 

read (22 ) ildim, j ldim, kldim 
read ( 2 2 ) duml , dum2 , dum3 , dum4 

read (22) ( ( (ql2 (i, 1, k, nx) , i=l, ildim) , k=l, kldim) , nx=l, 5) 

endif 

readj(5 ,*) dummy, nbstep 
if (iback .eq. 1) then 
back = .true, 
else 


c 

c find coordinates of center 
c 

xc = xstag + dels 
zc - zstag + dels * tan(alpr) 
yc - 0 . 
c 

sinthe = cos (alpr) 

costhe = sqrt (1 . -sinthe*sinthe) 

rr = dels/cosa 


c 

c starting location 
c 

c read(5, *,end=9999) isav, jsav, ksav, xf , yf , zf 

second = .false, 
igsav = 1 
issav = 1 
asav = 0. 
bsav = 0. 
gsav = 0 . 
ifatal = 0 
c 

c nsl is total number of streamlines to be generated 
c dtr is the conversion factor from degree to radian 
c dthe is the initial delta angle between two adjacent streamlines 
c 

dphi = angi 

if (nsl .gt. 1 .and. dphi . le . 0.) then 
dphi =180. / (nsl - 1) 
endif 

dphir = dphi * dtr 
cang = dphi * nsl 
jflag = 0 

if (cang .ge. 180.) jflag = 1 
endif 

c V 

c determin number of streamlines to be traced \ 0 

c 
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if(itype .ne. 0) then 

nsl * (ibmax-ibmin+1) * (kbmax-kbmin+1) 
istart * ibmin 
kstart * kbmin 
endif 

do 300 nl « l,nsl 
second = .false. 
if(iback .eq. 1) back = .true, 
if (back) then 

if(itype .eq. 0) then 
xin = xst 
yin = yst 
zin = zst 
else 

xin * xyz (istart/ 1, kstart, 1) 
yin * xyz {istart, 1, kstart, 2) 
zin * xyz (istart, 1, kstart, 3) 
endif 
else 

first « .false. 
if(nl .eq. 1) first = .true, 
c 

c initial position 
c 

if(nl .eq. 1) then 
xin = xstagl 
yin = ystagl 
zin — zstagl 

if (nsl .eq. 1 .and. jflag .eq. 1) then 
xin = xstag2 
yin * ystag2 
zin = zstag2 
endif 
else 

phir * (nl-l)*dphir 
cosp * cos (phir) 
sinp * sin (phir) 

write (6,*) 9 cosp, sinp cosp, sinp 
c 

c coordinate in (xc,yc,zc) system 
c 

xxp = -(rr * costhe) 
yyp = rr * sinthe * sinp 
zzp * -{rr * sinthe * cosp) 
write (6,*) 9 xxp,yyp,zzp 9 , xxp, yyp f zzp 
c 

c transfer coordinate to original (x,y, z) system 
c 

xin * xxp*cosa - zzp*sina + xc 
zin = xxp*sina + zzp*cosa + zc 
yin - yyp 

if (jflag .ne. 0 .and. nl .eq. nsl) then 
xin = xstag2 
yin * ystag2 
xin - zstag2 
if (alp .ge. 0.) then 
xin = xstagl 
yin = ystagl 
zin = zstagl 
endif 
endif 
endif 
endif 

c xin = xs (ispl+1, kdim-nl+1) 

? c* yin ^ys <ispl+l, kdim-nl+1) 

c k - *zin*= zs (ispl+1, kdim-nl+1) 



195 continue 

write (6,*) 9 xin,yin,zin 9 , xin, yin, zin 
ms * 1 
me = idim 
ns = 1 
ne * kdim 

call fndnrp (ms, me, ns,ne, xin, yin, zin) 
istat * 0 
delxi = 0 . 
delet = 0 . 

call quads (delxi , delet , xin, yin, zin, istat) 
if (istat .gt. 1) then 

write (6,*) ' Search failed, istat = 9 , istat 
write (40,*) 9 Location 9 , istart, kstart , 9 failed' 
go to 240 
endif 
il = isav 
kl = ksav 

c if (nl .eq. 1) then 

c if(kl .eq. 1 .or. kl .eq. kdim) then 

c delet = 0 . 

c endif 

c endif 

ell * delxi * delet 

clO = delxi - ell 

cOl = delet - ell 

cOO - 1. - delxi - delet + ell 

ipl = msp(il) 

kpl = nsp(kl) 

uin = u (il, 1, kl) *c00 + u (ipl, 1, kl) *cl0 + u (ipl, 1, kpl) *cll 

& + u (il, 1, kpl) *c01 

vin = v (il, 1, kl) *c00 + v (ipl, 1, kl) *cl0 + v (ipl, 1, kpl) *cll 

& + v (il, 1, kpl) *c01 

win = w (il, 1, kl) *c00 + w (ipl, 1, kl) *cl0 + w (ipl, 1, kpl) *cll 

& + w (il, 1, kpl) *c01 

pin = p (il, 1, kl) *c00 + p (ipl, 1, kl) *cl0 + p (ipl, 1, kpl) *cll 

& + p (il, 1, kpl) *c01 

dxdyin * dxdy (il, kl) *c00 + dxdy (ipl, kl) *clO + 

& dxdy (ipl, kpl) *cll + dxdy (il, kpl) *c01 

dxdzin = dxdz (il, kl) *c00 + dxdz (ipl, kl) *cl0 + 

& dxdz (ipl, kpl) *c 11 + dxdz (il, kpl) *c01 

c 

c find norm vector 


xnin = xn (il,kl) *c00 + xn (ipl, kl) *c!0 + 

& xn (ipl, kpl) *cll + xn (il, kpl) *c01 

ynin = yn (il, kl) *c00 + yn (ipl, kl) *cl0 + 

& yn (ipl, kpl) *cll + yn (il, kpl) *c01 

znin = zn (il,kl) *c00 + zn (ipl, kl) *cl0 + 

& zn (ipl, kpl) *cll + zn (il, kpl) *c01 

c 

c interpolate boundary layer transition parameters 

c 

if(flaura .eq. 9 y 9 .or. flaura .eq. 'Y') then 

repf = qll (il,l,kl,4)*c00 + qll (ipl, 1, kl, 4) *cl0 + 
& qll (ipl, l,kpl,4) *cll + qll (il, 1, kpl, 4) *c01 

c 

thkm = qll (il, l,kl,2) *c00 + qll (ipl, 1, kl, 2) *clO + 
& qll (ipl, 1, kpl, 2) *cll + qll (il, l,kpl, 2) *c01 

c 

thkd = qll (il, l,kl,3) *c00 + qll (ipl, 1, kl, 3) *cl0 + 
& qll (ipl, l,kpl, 3) *cll + qll (il, 1, kpl, 3) *c01 

c 

rek = ql2 (il,l,kl,3) *c00 + ql2 (ipl, 1, kl, 3) *cl0 + 

& q!2 (ipl, 1, kpl, 3) *cll + ql2 (il, 1, kpl, 3) *c01 



rem * repf * thkm / 12. 
red = repf * thkd / 12 . 
rex * repf * dist / 12. 
endif 
c 

x = 0 . 
y{l) = xin 
y (2) - yin 
y(3) * zin 

dist = sqrt ( (xin-xstag) **2+ (yin-ystag) **2+ (zin-zstag) **2) 
c 

c initial metric coefficient 
c 

hmet * dist 

c hmet * sqrt { (xin-xstag) **2+ (yin-ystag) **2+ (zin-zstag) **2) 

c 

c initial surface norm vector magnitude 
c 

surfn = sqrt(l. + dxdyin**2 + dxdzin**2) 
vel = sqrt (uin**2+vin**2+win**2 ) 
c 

c initial value for dydtau and dzdtau 
c 

y(4) = -hmet/surfn * (win/vel + uin/vel * dxdzin) 
y(5) *= hmet/surfn * (vin/vel + uin/vel * dxdyin) 
y (4) * -hmet * (xnin*win/vel - znin*uin/vel) 
y(5) * hmet * (xnin*vin/vel - ynin*uin/vel) 
c 

c modification of 8/16/93 
c 

y < 6 ) = hmet * (ynin*win/vel - znin*vin/vel) 
c end of modification 
c 

c re-compute hmet from dydtau and dzdtau 
c 

hmet = surfn/vel * (vin*y(5) - win*y(4)) 
hmet * (vin*y(5) -win*y (4) ) / (vel*xnin) 
c hmetl * -y(4) / (xnin*win/vel - znin*uin/vel) 

c hmet2 * y<5) / (xnin*vin/vel - ynin*uin/vel) 

hmet™ (y (5) -y (4) ) / (xnin* (vin+win) /vel-uin/vel* (ynin+znin) ) 
c 

c modification of 8/16/93 
c 

hmet ™ sqrt (y (4) **2 + y(5)**2 + y(6)**2) 
c end of modification 

c hmet = surfn/ (vel-uin/vel* (uin-vin*dxdyin-win*dxdzin) ) * 

c & (vin*y(5) - win*y(4) ) 

c hmetl = -surfn*y (4) / (win/vel+uin/vel*dxdzin) 

c hmet2 = surfn*y (5) / (vin/vel+uin/vel*dxdyin) 

c if(xf .ne. 0. .and. yf .ne. 0. .and. zf .ne. 0.) then 

C X = xf 

c y(l) * yf 

c y (2) - zf 

c endif 

c 

c output streamline locations 
c 

write (13,*) 9 Streamline No. ',nl 
write (14,*) ' Streamline No. ',nl 
write (6,*) 9 Streamline No. ',nl 
write (13, 9001) y (1) , y (2) ,y (3) , dist, hmet, pin 
if(flaura .eq. 9 y f .or. flaura .eq. 'Y') then 
write (26, 9001) thkm, thkd, rem, red, rek, rex 
endif 

write (14, 9002) dist/12 . ,hmet/12 . ,pin/pstag 

write (15)y(l),y(2),y(3) , -xnin, -ynin,-znin, vel, pin, hmet, isav,ksav 



c write (17,9010) y(l),y(2),y(3),y(l) -xnin,y (2) -ynin,y (3) -znin 

write (6, *) y(l),y(2),y(3),y(4),y(5) , surf n, vel, hmet, pin 
c 

c integration loop 
c 

h = hsave 
hmin * hmin 

write (6,*) 9 real step size', h, hmin 
xsave = y(l) 
ysave = y(2) 
zsave = y(3) 

c do 250 i = l,nstep 
200 continue 

if (second) then 
i = nstep 
else 

i = 1 
endif 

210 continue 

if (second) then 
y (1) - xb (i) 
y (2) - yb (i) 
y (3) - zb (i) 
if(i .eq. 2) then 
h = hb (i) 
else 

h - hb (i) - hb(i-l) 
endif 

c write(6,*) ' back ' , y (1) , y (2) , y (3) , h 

endif 

psave = surfp 

call rk34 (n, x, y, h, emax, hmin) 
if(ifatal .ne. 0) then 

write (40,*) ' Location ' , istart, kstart , ' failed' 
ifatal = 0 
go to 240 
endif 

if (back) then 

c if(itype .eq. 1) then 

distag = sqrt ( (y (1) -xstag) **2 + (y (2) -ystag) **2 
& + (y (3) -zstag) **2) 

if (distag .It. dels) then 
if (surfp .le. psave) then 
go to 220 

else if (abs (surfp-psave) .le. l.e-3) then 
go to 220 
endif 
endif 
i = i + 1 
nstep = i 
c else 

c if(i .gt. nstep) go to 220 

c i = i + 1 

c endif 

xb (i) = y (1) 
yb (i) = y (2 ) 
zb(i) = y ( 3) 
hb ( i) = x 
else 

i = i - 1 
endif 

c write (6,*) ' step ',i 

c write ( 6, *) isptr, xb (isptr) , yb (isptr) , zb (isptr) , hb (isptr) 

write (6,*) 'vp,y(5) ,wp,y (4) ' , vp,y (5) , wp,y (4) , surf n/ vel 
hmet = surfn / vel * (vp*y(5) - wp*y(4)) 
hmet = (vp*y(5) - wp*y (4) ) / ( vel*xnp) 



c 

c 


hmetl = -y (4) / (xnp*wp/vel - znp*up/vel) 
hmet2 = y(5) / (xnp*vp/vel - ynp*up/vel) 
hmet= (y (5) -y (4) ) / (xnp* (vp+wp) /vel-up/vel* (ynp+znp) ) 

c modification of 8/16/93 
c 

hmet = sqrt(y(4)**2 + y(5)**2 + y(6)**2) 
c end of modification 

c hmet = surfn/ (vel-up/vel* (up-vp*cLxdyp-wp*dxdzp) ) * 

c & (vp*y (5) - wp*y (4 ) ) 

c hmetl « -surfn*y (4) / (wp/vel+up/vel*dxdzp) 

c hmet2 = surfn*y (5) / (vp/vel+up/vel*dxdyp) 

dist - dist + sqrt { (y (1) -xsave) **2 + (y (2) -ysave) **2 
& + (y (3) -zsave) **2) 
xsave - y(l) 
ysave * y<2) 
zsave * y(3) 
hs = h 


c 

c 

c 


c 


c 


c 


c 


c 


interpolate boundary layer transition parameters 


if(flaura .eq. 
il = isav 
kl = ksav 


.or. flaura .eq. 9 Y' ) then 


ell = dxi * det 

clO = dxi - ell 

cOl « det - ell 

cOO * 1. - dxi - det + ell 

ipl * msp(il) 

kpl — nsp(kl) 

repf * qll (il,l,kl,4) *c00 + 
qll (ip lr If kpl, 4) *cll + 


qll (ipl, l,kl, 4)*cl0 + 
qll (il, l,kpl, 4) *c01 


thkm = qll (il,l,kl,2) *c00 + qll (ipl, 1, kl, 2) *clO + 
qll (ipl, 1, kpl, 2) *c 11 + qll (il, 1, kpl, 2) *c01 

thkd = qll (il,l,kl,3)*c00 + qll (ipl, l,kl, 3) *clO + 
qll (ipl, l,kpl, 3) *cll + qll (il, l,kpl, 3) *c01 

rek * ql2(il,l,kl,3)*c00 + ql2 (ipl, l,kl, 3) *cl0 + 
ql2 (ipl, 1, kpl, 3) *cll + ql2 (il, 1, kpl, 3) *c01 


rem = repf * thkm / 12 . 
red * repf * thkd / 12 . 
rex * repf * dist / 12 . 

write (26, 9001) thkm, thkd, rem, red, rek, repf 
endif 

write (13, 9001) y(l),y(2),y(3) , dist, hmet, surf p 
write (14, 9002) dist/ 12 ., hmet/ 12 . , surfp/pstag 

write (15) y(l),y(2),y{3) , -xnp, -ynp, -znp,vel, surf p, hmet, isav, ksav 
write (17, 9010) y(l),y(2),y(3) , y (1) -xnp, y (2) -ynp, y (3) -znp 
write (6,*) 9 i,k « ', isav, ksav 

write (6, *) y(l),y(2),y(3),y(4),y(5) , dist, hmet , surf p 
if ( . not . back ) then 

if(i .le. 1) go to 230 
endif 
go to 210 
continue 
if (back) then 
rewind 14 
rewind 13 
rewind 15 

if (flaura .eq. 9 y 9 .or. flaura .eq. 9 Y' ) then 
rewind 26 
endif 

dist=sqrt ( (y (1) -xstag) **2+ (y (2) -ystag) **2+ (y (3) -zstag) **2) 


220 


write (13,*) 9 Streamline No. ',nl 
write (14,*) r Streamline No. ',nl 
write (6,*) ' Streamline No. ',nl 
back = .false, 
second = .true, 
xin = xsave 
yin = ysave 
zin - zsave 
hsave - hs 
go to 195 
endif 

230 continue 
c 

c call subroutine heat to evaluate convective heating 
c 

call heat (qw, qwbtu, wltemp, cf ) 

rewind 13 

rewind 14 

rewind 15 

rewind 6 

rewind 24 

rewind 25 

write (30,*) istart, kstart, qw, qwbtu 
if (ifatal .ne. 0) then 

write (40,*) ' Location ', istart, kstart , 9 failed' 
ifatal = 0 
c go to 240 

endif 

if (itype .ne. 0) then 

q(istart, 1, kstart, 1) = qwbtu 
q ( istart, 1, kstart, 2) * wltemp 
q(istart, 1, kstart , 3) * deltal 
q(istart, 1, kstart, 4) = thetl 
q(istart, 1, kstart, 5) * cf 
endif 

240 continue 

if (istart .eq. ibmax) then 
if (kstart .eq. kbmax) then 
c 

c output plot3d file 
c 

jdim3d = 1 

if (.not. fext) then 

write (16) idim, jdim3d, kdim 

write (16) ( ( (xyz (i, j , k, 1) , i=l, idim) , j=l, jdim3d) , k=l, kdim) , 

& ( ( (xyz (i, j , k, 2 ) , i=l, idim) , j=l, jdim3d) , k=l, kdim) , 

& ( ( (xyz (i, j,k, 3) , i=l, idim) , j=l, jdim3d) ,k=l,kdim) 

endif 
c 

rewind 17 

write (17) idim, jdim3d, kdim 
write (17) aminf , alp, re, time 

write (17) ( ( ( (q(i, j,k, nx) , i*l, idim) , j=l, jdim3d) , k=l, kdim) , 

& nx=l,5) 

go to 9999 
else 

kstart = kstart + 1 
endif 
else 

kstart = kstart + 1 
if (kstart .gt. kbmax) then 
kstart = kbmin 
istart = istart + 1 
endif 
endif 

hsave = saveh 



300 continue 

if(itype .eq. 0 .and. .not. fext) then 
jdim3d * 1 

write (16) idim, jdim3d, kdim 

write (16) ( ( (xyz (i,j,k,l), i-1, idim) , j«l, jdim3d) , k*l, kdim) , 

& ( ( (xyz (i, j , k, 2 ) , i-1, idim) , j-1, jdim3d) ,k*l,kdim) , 

& ( ( (xyz (i,j,k,3),i=l, idim) , j=l, jdim3d) , k=l,kdim) 

c 

write ( 17 ) idim, jdim3d, kdim 
write (17) aminf , alp, re , time 

write (17) ( ( ( (q(i, j,k,nx) , i-1, idim) , j=l, jdim3d) ,k=l,kdim) , 

& nx=l,5) 

endif 
c 

9001 format (6 (lx, ell . 4, lx) ) 

9002 format (3 (lx, el5. 6, lx)) 

9010 format (6 (el2 .5, lx) ) 

9999 continue 
stop 
end 
c 

c subroutine rk34 
c 

subroutine rk34 (n, x, y, h, e, hmin) 
c 

common /save/ igsav, issav, isav, jsav,ksav, asav,bsav,gsav, 

& first, back, second, ifatal 

c 

dimension y(50), ytemp(50), yhat(50), f(50), a (50), 

&b (50) , c (50) , d (50) 
logical second, first , back 
c 

data cO, c2, c3, chO, ch2, ch3, ch4, al, a2, a3,bl0,b20,b21, 

& b30,b31,b32,b40,b42,b43/. 16122448 98, .5998345284, 

& .2389409818, .1557823129, .6205184777, .1681436539, 

& .0555555556, .2857142857, .4666666667, .9210526316, 

& .2857142857, .0855555556, .3811111111, 

& .5574792244, -1.406455023, 1.770028430, .1612244898, 

& .5998345284, .2389409818/ 

c 

1 continue 
xtemp = x 

do 2 i s l,n 
ytemp(i) - y(i) 

2 continue 

call funct (xtemp, ytemp, f) 

if (ifatal .ne. 0) return 

xtemp = x + al*h 

do 3 i - l,n 

a(i) * h * f(i) 

ytemp(i) - y(i) + bl0*a(i) 

3 continue 

call funct (xtemp, ytemp, f ) 
if (ifatal .ne. 0) return 
xtemp = x + a2*h 
do 4 i = l,n 
b (i) - h * f (i) 

ytemp(i) * y(i) + b20*a(i) + b21*b(i) 

4 continue 

call funct (xtemp, ytemp, f) 
if (ifatal .ne. 0) return 
xtemp = x + a3*h 
do 5 i * l,n 
c (i) = h * f (i) 

ytemp(i) « y(i) + b30*a(i) + b31*b(i) + b32*c(i) 

5 continue 



6 


call funct (xtemp, ytemp, f ) 
if {if atal .ne. 0) return 
xtemp - x + h 
do 6 i = l,n 
d(i) = h * f (i) 

ytemp(i) - y(i) + b40*a(i) + b42*c(i) + b43*d<i) 
continue 

call funct (xtemp , ytemp, f) 
if{ifatal .ne. 0) return 
do 7 i = l,n 

ytemp(i) = y(i) + c0*a(i) + c2*c(i) + c3*d(i) 
yhat(i) = y(i) + ch0*a(i) + ch2*c(i) + ch3*d(i) + ch4*h*f(i) 
a(i) = abs(yhat(i) - ytemp (i) ) 

7 continue 

do 8 i * l,n 
c(i) = abs(yhat(i)) 
if (c (i) .le. a(i)*l.e-3) c(i) = 1. 
c if(abs(f(i)) .le. l.e-5) c(i) = 1. 

b(i) = a (i) / c ( i ) 

if (.not. second .and. b(i) .gt. e) go to 11 

8 continue 
c 

9 continue 
x = x + h 
etemp = e / 16. 
iflag = 0 

do 10 i = l,n 
y (i) = yhat (i) 

if(b(i) .gt. etenp) iflag = 1 

10 continue 

call funct (x,y r f) 
if(ifatal .ne. 0) return 
if (iflag .eq. 1) return 
h = h + h 

c if(h .gt. 1.0e~3) h = 1.0e-3 

return 

11 continue 

if(abs(h) .gt. abs (hmin) ) go to 12 

write (6,*) 9 relative truncation error criterion could 
& not be satisfied/' , e, hmin 
c stop 

go to 9 

12 continue 

h = h * 0.5 
write (6 , *) 9 h = ' , h 
go to 1 
end 


c 

c subroutine funct 


c 


c 

c 

c 


& 

& 

& 

& 


subroutine funct (x,y,f) 

common /nomxyz/ xn(151,110),yn(151,110),zn (151,110)/ xnp, ynp, znp 
common /sfbond / msp (151) ,msr (151) , nsp (151) ,nsr (151) 
common /grdxyz/ xyz (151/ 40, 110, 3) , idim, jdim, kdim, 

iblank(151, 40, 110) , blank, isubs (2 ) , 
jsubs (2) , ksubs (2 ) 

common /qxyz/ u (151, 40, 110) , v ( 151, 40, 110) , 
w (151, 40, 110) , p (151, 40, 110) 

common /dvdwyz/ dvdy (151, 110) ,dvdz (151, 110) ,dwdy (151, 110) , 
dwdz (151, 110) 


modification of 8/16/93 


common /duxyz/ dudx (151, 110) ,dudy (151,110) ,dudz (151,110) , 
& dvdx (151,110) ,dwdx(151, 110) 

end of modification 


c 



common /geom/ dxdy (151,110) , dxdz (151, 110) 

common / surfpr/ vel,surfn,up,vp,wp,surfp,dxdzp,dxdyp,dxi,det 
common / save/ igsav, issav, isav, jsav,ksav,asav,bsav,gsav, 

& first, back, second, ifatal 

logical first, back, second 
c 

dimension y(50),f(50) 
c 

sx * y(l) 
sy « y (2) 
sz « y(3) 
c 

c surface interpolation 
c 

ms = isav - 3 
me = isav + 3 
if (ms .It. 1) ms = 1 
if (me . gt . idim) me = idim 
ns = ksav - 3 
ne = ksav + 3 
if (ns .It. 1) ns = 1 
if (ne .gt. kdim) ne = kdim 
c ns *= 1 

c ne = kdim 

call fndnrp (ms, me, ns, ne, sx, sy, sz) 
c 

istat * 0 
delxi * 0 . 
delet = 0 . 

call quads (delxi, delet , sx, sy, sz, istat) 
if (istat .gt. 1) then 

write (6,*) 9 Search failed, istat = istat 
ifatal * 1 
return 
endif 
il = isav 
kl = ksav 

c if (first) then 

c if (kl .eq. 1 .or. kl .eq. kdim) then 

c delet = 0 . 

c endif 

c endif 

dxi = delxi 

det * delet 

ell * delxi * delet 

clO * delxi - ell 

cOl = delet - ell 

cOO * 1. - delxi - delet + ell 

ipl = msp(il) 

kpl = nsp(kl) 

up - u (il, 1, kl) *c00 + u (ipl, l,kl) *cl0 + u (ipl, 1, kpl) *cll 
& + u (il, 1, kpl) *c01 

vp - v (il, 1, kl) *c00 + v (ipl, l,kl) *cl0 + v (ipl, 1, kpl) *cll 
& + v (il, 1, kpl) *c01 

wp = w(il, l,kl) *c00 + w (ipl, l,kl) *cl0 + w (ipl, 1, kpl) *cll 
& + w(il, l,kpl) *c01 

surfp = p (il, l,kl) *c00 + p (ipl, 1, kl) *cl0 + p (ipl, 1, kpl) *cll 
& + p (il, 1, kpl) *c01 

dvdyp = dvdy (il, kl) *cQ0 + dvdy (ipl, kl) *cl0 + 

& dvdy (ipl, kpl) *cll + dvdy (il, kpl) *c01 

dvdzp = dvdz (il, kl) *c00 + dvdz (ipl, kl) *cl0 + 

& dvdz (ipl, kpl ) *cll + dvdz (il, kpl) *c01 

dwdyp = dwdy (il, kl) *c00 + dwdy (ipl, kl) *cl0 + 

& dwdy (ipl, kpl) *cll + dwdy (il, kpl) *c01 

dwdzp = dwdz (il, kl) *c00 + dwdz (ipl, kl) *cl0 + 

& dwdz (ipl, kpl) *cll + dwdz (il, kpl) *c01 



c 

c modification of 8/16/93 
c 

dudxp = dudx (il, kl) *c00 + dudx (ipl,kl) *clO + 

& dudx(ipl, kpl) *cll + dudx (il, kpl) *c01 

dudyp = dudy (il, kl) *c00 + dudy (ipl, kl) *clO + 

& dudy (ipl, kpl) *cll + dudy (il, kpl) *c01 

dudzp * dudz (il, kl) *c00 + dudz (ipl,kl) *clO + 

& dudz (ipl, kpl) *cll + dudz (il, kpl) *c01 

dvdxp * dvdx (il, kl) *c00 + dvdx (ipl, kl) *clO + 

& dvdx (ipl, kpl) *cll + dvdx (il, kpl) *c01 

dwdxp = dwdx (il, kl) *c00 + dwdx (ipl, kl) *clO + 

& dwdx (ipl, kpl) *cll + dwdx (il, kpl) *c01 

c 

dxdyp = dxdy (il,kl) *c00 + dxdy (ipl, kl) *clO + 

& dxdy (ipl, kpl) *cll + dxdy (il, kpl) *c01 

dxdzp = dxdz (il,kl) *c00 + dxdz (ipl, kl) *clO + 

& dxdz (ipl, kpl) *cll + dxdz (il, kpl) *c01 

c 

c surface norm vector 
c 

xnp = xn (il, kl) *c00 + xn (ipl, kl) *clO + 

& xn (ipl, kpl) *cll + xn (il, kpl) *c01 

ynp = yn (il, kl) *c00 + yn (ipl, kl) *clO + 

& yn (ipl, kpl ) *cll + yn (il, kpl) *c01 

znp « zn (il, kl) *c00 + zn (ipl, kl) *clO + 

& zn(ipl,kpl) *cll + zn (il,kpl) *c01 

c 

vel = sqrt (up*up + vp*vp + wp*wp) 
surfn = sqrt(l. + dxdyp * dxdyp + dxdzp*dxdzp) 
c 

f ( 1 ) * up 
f (2) = vp 
f (3) = wp 

c f ( 4 ) = dvdyp * y ( 4 ) + dvdzp * y ( 5 ) 

c f(5) = dwdyp * y(4) + dwdzp * y{5) 

f ( 4 ) = dvdxp * y ( 6 ) + dvdyp * y ( 4 ) + dvdzp * y ( 5 ) 

f(5) = dwdxp * y(6) + dwdyp * y ( 4 ) + dwdzp * y(5) 

f(6) * dudxp * y(6) + dudyp * y(4) + dudzp * y(5) 

if (back) then 
f(l) = -f(l) 
f (2 ) = -f (2 ) 
f (3) = -f (3) 
endif 

c write ( 6, *) r f(l),f(2) - ' , f ( 1) , f (2) , f (3) , f (4) , f (5) 

return 
end 
c 

subroutine bounds ( mdim, ndim,ms,me, ns, ne, 

& msper,nsper) 

common /sfbond/ msp (151) ,msr (151) ,nsp (151) ,nsr (151) 
c 

c Initialize arrays, 
c 

do 10 m = 1 , mdim 
msp (m) = 0 

msr(m) = 0 

10 continue 

do 11 n = l,ndim 
nsp(n) = 0 

nsr(n) = 0 

11 continue 
c 

c Set up +/- index arrays, 
c 

do 20 m = ms, me 



20 


msp (m) * m+1 
msr(m) - m-1 
continue 

if (msper.eq.O) then 
msp (me) = me 
msr (ms) = ms 
else 

msp (me) = ms+1 
msr (ms) = me-1 
endif 

do 21 n = ns,ne 
nsp(n) = n+1 
nsr(n) * n-1 
21 continue 

if (nsper.eq.O) then 
nsp(ne) = ne 
nsr(ns) = ns 
else 

nsp(ne) = ns+1 
nsr(ns) = ne-1 
endif 

return 

end 

c 

subroutine surmet ( mdim, ndim,ms, me, ns,ne,msper, nsper) 
c 

c Compute metrics for defining surface grid, 
c 

common /nomxyz/ xn (151, 110) , yn (151, 110) , zn (151# 110) ,xnp, ynp, znp 
common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151, 40,110), blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

common /sfxiet/ xixs (151, 110) ,xiys (151, 110) ,xizs (151, 110) , 

& etxs (151,110) , etys (151, 110) ,etzs (151, 110) , 

& flgsm(151) , f lgsn (151) 

common /sfbond/ msp{151) ,msr (151) ,nsp(151) ,nsr (151) 
common /qxyz/ u (151, 40, 110) , v (151, 40, 110) , 

& w(151, 40,110) , p(151, 40,110) 

common /dvdwyz/ dvdy (151, 110) ,dvdz (151, 110) ,dwdy (151, 110) , 

& dwdz (151, 110) 

c 

c modification of 8/16/93 

c 

common /duxyz/ dudx(151, 110) ,dudy (151, 110) ,dudz (151, 110) , 

& dvdx (151, 110) , dwdx (151, 110) 

c end of modification 

common /geom/ dxdy (151, 110) ,dxdz (151, 110) 
dimension jacm(151) , jacn (151) 
c 

njac = 0 
c 

c Initialize arrays, 
c 

do 10 m = l,mdim 
flgsm(m) = 0. 
jacm(m) = 0 
jacn(n) = 0 

10 continue 

do 11 n = l,ndim 
flgsn(n) = 0. 

11 continue 

do 12 m = l,mdim 
do 12 n = l,ndim 
xixs (m,n) = 0. 



xiys(m, n) = 0. 
xizs(m, n) = 0. 
etxs(m,n) = 0. 
etys(m, n) = 0. 
etzs (m, n) = 0 . 
dvdy(m,n) = 0. 
dvdz(m, n) = 0. 
dwdy(m, n) = 0. 
dwdz(m,n) = 0. 
dxdy (m, n) = 0. 
dxdz(m, n) = 0. 

12 continue 
c 

c Set up weight functions, 
c 

do 20 m * ms , me 
f lgsm <m) = .5 

20 continue 

if (msper.eq.O) then 
flgsm(me) = 1. 
f lgsm (ms) - 1. 

endif 

do 21 n = ns,ne 
flgsn(n) = .5 

21 continue 

if (nsper.eq.0) then 
flgsn(ns) = 1. 
flgsn(ne) = 1. 

endif 

c 

c Compute metrics, 
c 

do 30 m * ms , me 
mp = msp (m) 
mr *= msr (m) 
do 30 n = ns,ne 
np *= nsp(n) 
nr = nsr(n) 

xxi - (xyz (mp, l,n, 1) - xyz (mr, 1, n, 1) ) * flgsm(m) 

yxi = (xyz (mp, 1, n, 2) - xyz (mr, l,n, 2) ) * flgsm(m) 

zxi = (xyz (mp, 1, n, 3) - xyz (mr, 1, n, 3) ) * flgsm(m) 

vxi = (v(mp,l,n) - v(mr,l,n)) * flgsm(m) 

wxi = (w(mp,l,n) - w(mr,l,n)) * flgsm(m) 

xet = (xyz (m, l,np, 1) - xyz (m, 1, nr, 1) ) * flgsn(n) 

yet = (xyz (m, l,np, 2) - xyz (m, 1, nr, 2) ) * flgsn(n) 

zet * (xyz (m, l,np, 3) - xyz (m, 1, nr, 3) ) * flgsn(n) 

vet - (v(m, l,np) - v (m, l,nr) ) * flgsn(n) 

wet = (w(m, l,np) - w(m, l,nr) ) * flgsn(n) 

c 

c modification of 8/16/93 

c 

uxi = (u(mp,l,n) - u(mr,l,n)) * flgsm(m) 
uet = (u(m, l,np) - u (m, l,nr)) * flgsn(n) 
c end of modification 

c 

c debug print 
c 

if (n .eq. ns .or. n .eq. ne) then 
write (48, *) m, n, xxi, yxi, zxi 
write (49, *) m, n, xet, yet, zet 
endif 
c 

c confute surface norm 
c 


if (m .eq. ms) then 



xn (m, n) = 1. 
yn(m, n) * 0. 
zn(m, n) = 0. 

else if (m .eq. me .and. n .ne. ns .and. n ,ne. ne)then 
xl = xyz (m-1, l,n+l, 1) - xyz (m, l,n,l) 
yl = xyz (m-1, 1, n+1, 2) - xyz(m,l,n,2) 
zl = xyz (m-1, l,n+l, 3) - xyz (m, l,n,3) 
x2 = xyz (m, 1, n+1 , 1 ) - xyz (m-1, 1, n, 1) 
y2 = xyz (m, 1, n+1, 2) - xyz (m-1, 1, n, 2) 
z2 * xyz (m, 1, n+1 , 3 ) - xyz (m-1, 1, n, 3) 
xnl = yl*z2 - y2*zl 
ynl = x2*zl - xl*z2 
znl = xl*y2 - x2*yl 

sqri *= 1 . /sqrt (xnl*xnl + ynl*ynl + znl*znl) 
xnl = xnl*sqri 
ynl ■= ynl*sqri 
znl * znl*sqri 

xl = xyz (m-1, 1, n-1, 1) - xyz(m, l,n,l) 
yl * xyz (m-1, 1, n-1, 2) - xyz(m, l,n,2) 
zl = xyz (m-1, l,n-l, 3) - xyz (m, l,n,3) 
x2 = xyz (m-1, l,n, 1) - xyz (m, 1, n-1, 1) 
y2 = xyz (m-1, l,n,2) - xyz (m, 1, n-1, 2) 
z2 = xyz (m-1, 1, n, 3) - xyz (m, 1, n-1, 3) 
xn2 = yl*z2 - y2*zl 
yn2 - x2*zl - xl*z2 
zn2 = xl*y2 - x2*yl 

sqri = 1 . /sqrt (xn2*xn2 + yn2*yn2 + zn2*zn2) 
xn2 =* xn2*sqri 
yn2 = yn2*sqri 
zn2 = zn2*sqri 

xnn = 0 .5* (xnl+xn2) 
ynn * 0 . 5* (ynl+yn2 ) 
znn = 0 . 5* (znl+zn2) 

sqri = 1 . /sqrt (xnn *xnn + ynn*ynn + znn*znn) 
xn(m,n) = xnn*sqri 
yn(m,n) = ynn*sqri 
zn(m, n) = znn*sqri 

else if (n .eq. ns .and. m .ne. ms .and. m .ne. me) then 
xl = xyz (m-1, 1, n+1, 1) - xyz (m, l,n,l) 
yl = xyz (m-1, 1, n+1, 2) - xyz(m,l,n,2) 
zl = xyz (m-1, 1, n+1, 3) - xyz(m,l,n,3) 
x2 » xyz (m, l,n+l, 1) - xyz (m-1, 1, n, 1) 
y2 = xyz (m, 1, n+1, 2) - xyz (m-1, 1, n, 2) 
z2 = xyz (m, 1, n+1, 3) - xyz (m-1, l,n, 3) 
xnl = yl*z2 - y2*zl 
ynl = x2*zl - xl*z2 
znl — xl*y2 - x2*yl 

sqri = 1 . /sqrt (xnl *xnl + ynl*ynl + znl*znl) 
xnl = xnl*sqri 
ynl = ynl*sqri 
znl = znl*sqri 

xl = xyz (m+1, 1, n+1, 1) - xyz (m, l,n,l) 
yl = xyz (m+1, 1, n+1, 2) - xyz(m,l,n,2) 
zl = xyz (m+1, 1, n+1, 3) - xyz (m, l,n,3) 
x2 = xyz (m+1, l,n, 1) - xyz (m, 1, n+1 , 1) 
y2 * xyz (m+1, l,n, 2) - xyz (m, 1, n+1, 2) 
z2 = xyz (m+1, l,n, 3) - xyz (m, 1, n+1, 3) 
xn2 = yl*z2 - y2*zl 
yn2 = x2*zl - xl*z2 
zn2 = xl*y2 - x2*yl 

sqri « 1. /sqrt (xn2*xn2 + yn2*yn2 + zn2*zn2) 
xn2 = xn2*sqri 
yn2 = yn2*sqri 



zn2 = zn2*sqri 
xnn = 0 . 5* (xnl+xn2) 
ynn = 0 . 5* (ynl+yn2) 
ynn * 0 . 

znn - 0 . 5* (znl+zn2) 

sqri * 1 . /sqrt (xnn*xnn + ynn*ynn + znn*znn) 
xn (m, n) * xnn*sqri 
yn(m,n) = ynn*sqri 
zn(m, n) = znn*sqri 

else if(n .eq. ne .and. m .ne. ms .and. m .ne. me) then 
xl = xyz (m+1, l,n-l, 1) - xyz(m, l,n,l) 
yl = xyz (m+1, 1, n-1, 2) - xyz (m, l,n,2) 
zl = xyz (m+1, l,n-l, 3) - xyz(m, l,n,3) 
x2 = xyz (m, 1, n-1, 1) - xyz (m+1, 1, n, 1) 
y2 - xyz (m, 1, n-1, 2) - xyz (m+1, 1, n, 2 ) 
z2 = xyz (m, 1, n-1, 3) - xyz (m+1, 1, n, 3) 
xnl = yl*z2 - y2*zl 
ynl = x2*zl - xl*z2 
znl = xl*y2 - x2*yl 

sqri = 1 . /sqrt (xnl*xnl + ynl*ynl + znl*znl) 
xnl = xnl * sqri 
ynl = ynl *sqri 
znl = znl*sqri 

xl = xyz (m-1, 1, n-1, 1) - xyz(m, l,n,l) 
yl = xyz (m-1, 1, n-1, 2) - xyz(m r l/n/2) 
zl = xyz (m-1/ 1/ n-1/ 3) - xyz(m/l/n,3) 
x2 * xyz (m-1/ 1/ n/ 1) - xyz (m r 1/ n-1/ 1) 
y2 * xyz (m-1, 1, n, 2) - xyz (m f 1, n-1/ 2) 
z2 * xyz (m-1/ 1, n, 3) - xyz (m, 1/ n-1, 3) 
xn2 = yl*z2 - y2*zl 
yn2 = x2*zl - xl*z2 
zn2 = xl*y2 - x2*yl 

sqri = 1 , /sqrt (xn2*xn2 + yn2*yn2 + zn2*zn2) 

xn2 = xn2*sqri 

yn2 = yn2*sqri 

zn2 = zn2*sqri 

xnn » 0 . 5* (xnl+xn2) 

ynn = 0 . 5* (ynl+yn2 ) 

ynn - 0 . 

znn = 0 . 5* <znl+zn2) 

sqri = 1 . /sqrt (xnn* xnn + ynn*ynn + znn*znn) 
xn(m,n) = xnn*sqri 
yn(rri/n) = ynn*sqri 
zn(m,n) * znn*sqri 

else if (m .eq. me .and. n. eq. ns) then 
xl = xyz (m-1/ 1/ n, 1) - xyz(m r l/n/l) 
yl * xyz (m-1/ 1/ n, 2) - xyz(m,l,n,2) 
zl = xyz (m-1/ l/n f 3) - xyz(m,l,n,3) 
x2 = xyz (m, l,n+l, 1) - xyz(m/l/n/l) 
y2 = xyz (m, l,n+l, 2) - xyz (m, 1, n, 2) 
z2 * xyz (m, l,n+l, 3) - xyz(m/l f n,3) 
xnl = yl*z2 - y2*zl 
ynl = x2*zl - xl*z2 
ynl - 0 . 

znl = xl*y2 - x2*yl 

sqri = 1 . /sqrt (xnl *xnl + ynl*ynl + znl*znl) 
xn(m,n) = xnl*sqri 
yn (m, n) «= ynl*sqri 
zn(m,n) = znl*sqri 

else if (m .eq. me .and. n .eq. ne) then 
xl = xyz (m f l,n-l/ 1) - xyz(ni/l/n/l) 
yl = xyz (m f 1, n-1, 2) - xyz(m,l,n,2) 
zl = xyz (m, 1, n-1, 3) - xyz (m, 1, n, 3) 
x2 = xyz (m-1, 1, n, 1) - xyz(m/l f n/l) 
y2 = xyz (m-1, 1, n, 2) - xyz(m,l,n/2) 



z2 = xyz (m-1, 1, n, 3) - xyz (m, l,n,3) 
xnl = yl*z2 - y2*zl 
ynl * x2*zl - xl*z2 
ynl = 0 . 

znl = xl*y2 - x2*yl 

sqri = 1 . /sqrt (xnl*xnl + ynl*ynl + znl*znl) 
xn(m,n) = xnl*sqri 
yn (m, n) = ynl*sqri 
zn(m, n) = znl*sqri 
else 

xl * xyz (m-1, l,n-l, 1) - xyz (ro, 1, n, 1) 
yl * xyz (m-1, l,n-l, 2) - xyz (m, 1 , 11 , 2 ) 
zl = xyz (m-1, l,n-l, 3) - xyz (m, l,n,3) 
x2 = xyz (m-1, l,n, 1) - xyz (m, 1, n-1, 1) 
y2 = xyz (m-1, l,n,2) - xyz (m, 1, n-1, 2) 
z2 = xyz (m-1, l,n, 3) - xyz (m, 1, n-1, 3) 
xnl * yl*z2 - y2*zl 
ynl = x2*zl - xl*z2 
znl = xl*y2 - x2*yl 

sqri *= 1 . /sqrt (xnl*xnl + ynl*ynl + znl*znl) 
xnl - xnl*sqri 
ynl = ynl*sqri 
znl = znl*sqri 

xl * xyz (m-1, 1, n+1, 1) - xyz (m, l,n,l) 
yl = xyz (m-1, 1, n+1, 2) - xyz (m, l,n,2) 
zl = xyz (m-1, 1, n+1, 3) - xyz (m, l,n,3) 
x2 = xyz (m, l,n+l, 1) - xyz (m-1, 1, n, 1) 
y2 * xyz (m, 1, n+1, 2) - xyz (m-1, 1, n, 2) 
z2 * xyz (m, 1, n+1, 3) - xyz (m-1, 1, n, 3) 
xn2 = yl*z2 - y2*zl 
yn2 = x2*zl - xl*z2 
zn2 = xl*y2 - x2*yl 

sqri = 1 , /sqrt (xn2*xn2 + yn2*yn2 + zn2*zn2) 
xn2 = xn2 *sqri 
yn2 = yn2*sqri 
zn2 = zn2*sqri 

xl = xyz (m+1, 1, n+1, 1) - xyz (m, l,n,l) 
yl = xyz (m+1, 1, n+1, 2) - xyz (m, l,n,2) 
zl — xyz (m+1, 1, n+1, 3) - xyz (m, l,n,3) 
x2 - xyz (m+1, 1, n, 1) - xyz (m, 1, n+1, 1) 
y2 - xyz (m+1, l,n,2) - xyz (m, 1, n+1, 2) 
z2 - xyz (m+1, l,n, 3) - xyz (m, 1, n+1, 3) 
xn3 = yl*z2 - y2*zl 
yn3 = x2*zl - xl*z2 
zn3 * xl*y2 - x2*yl 

sqri = 1 . /sqrt (xn3*xn3 + yn3*yn3 + zn3*zn3) 
xn3 = xn3*sqri 
yn3 = yn3*sqri 
zn3 = zn3*sqri 

xl = xyz (m+1, 1, n-1, 1) - xyz (m, l,n,l) 
yl = xyz (m+1, 1, n-1, 2) - xyz(m,l,n,2) 
zl = xyz (m+1, 1, n-1, 3) - xyz (m, l,n,3) 
x2 = xyz (m, 1, n-1, 1) - xyz (m+1, 1, n, 1) 
y2 = xyz (m, 1, n-1, 2) - xyz (m+1, 1, n, 2) 
z2 *= xyz (m, 1, n-1, 3) - xyz (m+1, 1, n, 3) 
xn4 = yl*z2 - y2*zl 
yn4 = x2*zl - xl*z2 
zn4 = xl*y2 - x2*yl 

sqri = 1 . /sqrt (xn4*xn4 + yn4*yn4 + zn4*zn4) 
xn4 = xn4*sqri 
yn4 = yn4*sqri 
zn4 = zn4*sqri 



xnn = 0 .25* (xnl+xn2+xn3+xn4) 
ynn * 0 .25* (ynl+yn2+yn3+yn4) 
znn = 0 . 25* (znl+zn2+zn3+zn4) 
sqri = 1 . /sqrt (xnn*xnn + ynn*ynn + znn*znn) 
xn(m, n) = xnn* sqri 
yn(m, n) = ynn* sqri 
zn(m, n) = znn*sqri 
endif 
c 

rxrx = xxi**2 + yxi**2 + zxi**2 

rxre = xxi*xet + yxi*yet + zxi*zet 

rere - xet**2 + yet**2 + zet**2 

rjac * ( rxrx*rere - rxre**2 ) 
c 

c Make sure jacobian is positive. For degenerate or collapsed points 
c a more rigorous forward or backward scheme should be used. For now, 
c just ignore these points, 
c 

if (rjac.gt.O. .and. rjac .gt. l.e-10) then 
c if (rjac .gt. 0.) then 

rd - 1 . / rjac 

xixs(m,n) = ( rere*xxi - rxre*xet) * rd 

xiys(m,n) = ( rere*yxi - rxre*yet) * rd 

xizs(m,n) = ( rere*zxi - rxre*zet) * rd 

etxs (m,n) = (-rxre*xxi + rxrx*xet) * rd 

etys(m, n) = <~rxre*yxi + rxrx*yet) * rd 

etzs (m,n) = (-rxre*zxi + rxrx*zet) * rd 

c 

c surface velocity derivatives 
c 

dvdy (m, n) = vxi*xiys (m, n) + vet*etys (m, n) 
dvdz (m, n) = vxi*xizs (m, n) + vet*etzs (m, n) 
dwdy(m, n) = wxi*xiys (m, n) + wet*etys (m, n) 
dwdz (m, n) = wxi*xizs (m, n) + wet*etzs (m, n) 
c 

c modification of 8/16/93 
c 

dudx(m, n) * uxi*xixs (m, n) + uet*etxs <m f n) 
dudy(m, n) = uxi*xiys (m f n) + uet*etys (m, n) 
dudz (m, n) = uxi*xizs (m f n) + uet*etzs (m, n) 
dvdx(m f n) = vxi*xixs (m, n) + vet*etxs (m r n) 
dwdx (m, n) = wxi*xixs (m r n) + wet*etzs (m r n) 
c end of modification 
c 

c surface derivatives 
c 

dxdy(m, n) * xxi*xiys (m f n) + xet*etys (m, n) 
dxdz(m, n) = xxi*xizs (m, n) + xet*etzs (m, n) 
c if (n .eq. ns .or. n .eq. ne) then 

c write (50, *) m,n,dvdy (m,n) , dvdz (m,n) 

c write (51,*) m, n, dwdy (m, n) , dwdz (m, n) 

c write (52, *) m, n, dxdy (m, n) , dxdz (m, n) 

c endif 

else 

njac = njac + 1 
jacm(njac) = m 
jacn(njac) - n 
xixs (m, n) = 0 . 
xiys (m, n) = 0 . 
xizs (m, n) = 0 . 
etxs (m, n) = 0 . 
etys (m, n) = 0. 
etzs (m,n) = 0. 
dvdy(m, n) = 0. 
dvdz (m, n) = 0. 
dwdy (m, n) * 0. 
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dwdz (m, n) - 0. 
dxdy(m, n) * 0. 
dxdz (m, n) = 0. 
endif 

continue 
c 

c Write questionable grid points to output log file, 
c 

if (njac.ne.O) then 

open (2, f ile=' surmet . log' , status*' unknown' , form=' formatted' ) 
do 40 n = l,njac 

write (2,1010) jacm(n) , jacn (n) 

40 continue 

close (2) 
endif 

1010 format ('Metric computation for grid point (' , i2, ' , ' , i2, 
failed. This point will be ignored.') 
return 
end 
c 

subroutine fndnrp { ms, me, ns, ne, 

& xp,yp,zp) 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim,kdim, 

& iblank (151, 40, 110) , blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

common /save/ igsav, issav,msave, jsav,nsave,asav,bsav,gsav, 

& first, back, second, ifatal 

logical blank, first, back, second 
c 

c Find point on surface subset that is closest to (xp,yp,zp). 
c 

dsave ■= 1.0e8 
do 10 n = ns,ne 

do 10 m * ms, me 

dist = (xyz (m, l,n, 1) -xp) **2 + 

1 (xyz (m, l,n,2) -yp) **2 + 

2 (xyz (m, 1, n, 3) -zp) **2 
if (dist . It .dsave) then 

dsave = dist 
msave = m 
nsave = n 
endif 

10 continue 

c write (6,*) 9 msave, nsave * ', msave, nsave 

return 
end 
c 

subroutine quads ( delxi, delet , xp, yp, zp, istat ) 
c 

parameter (toler=0 . 1) 

common /save/ igsav, issav, js, jsav,ks,asav,bsav,gsav, 

& first, back, second, ifatal 

common / grdxyz/ xyz (151, 40, 110, 3) , jsmax, ismax, ksmax, 

& iblank (151,40,110), blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

common /sfxiet/ xixs (151,110) , xiys (151,110), xizs (151,110), 

& etxs (151, 110) ,etys (151, 110) ,etzs (151, 110) , 

& flgsm(151) ,flgsn(151) 

common /sfbond/ jsp(151) ,jsr(151) ,ksp(151) , ksr (151) 
logical first, back, second 
ifail = 0 
c 

c find quadrant in which the point (xp,yp,zp) lies, 
c 



jstart = js 
kstart * ks 

5 continue 

iterm = ((jsmax + ksmax)/2) + 2 
do 10 iter = 1, iterm 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 
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c 

c 

c 
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& 

& 

& 

& 

1 

1 

1 

1 

1 
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1 

1 

1 

1 

1 

1 
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delxi = 0.5 
delet * 0.5 
if(js .eq. 1) then 
delxi = .25 
delet = .25 
endif 

do 11 it *= 1,10 
if(js .eq. 1) then 

ell = 1 - delxi - delet 
cl2 = delet 
cl3 = delxi 
jp = jsp(js) 
kp = ksp(ks) 

xt = xyz ( js, 1, ks, 1) *cll + xyz ( jp, 1 , ks, 1) *cl2 
+ xy z ( j p , 1 , kp , 1 ) * c 1 3 

yt *= xyz ( js, 1, ks, 2) *cll + xyz ( jp, 1, ks, 2) *cl2 
+ xyz ( jp, l,kp,2) *cl3 

zt = xyz ( js, l,ks, 3) *cll + xyz ( jp, 1, ks, 3) *cl2 
+ xyz ( jp, 1, kp, 3) *cl3 

xt = xyz ( js, 1, ks, 1) *cll + xyz ( jp, l,ks, 1) *cl2 
+ xyz ( jp, 1, kp, 1) *cl3 

xix = xixs(js,ks) * ell + xixs(jp,ks) * cl2 
+ xixs(jp,kp) * cl3 

xiy = xiys(js,ks) * ell + xiys(jp,ks) * cl2 
+ xiys<jp,kp) * cl3 

xiz = xizs(js,ks) * ell + xizs(jp,ks) * cl2 
+ xizs(jp,kp) * cl3 

etx = etxs(js,ks) * ell + etxs(jp,ks) * cl2 
+ etxs(jp,kp) * cl3 

ety = etys(js,ks) * ell + etys(jp,ks) * cl2 
+ etys ( jp, kp) * cl3 

etz = etzs(js,ks) * ell + etzs(jp,ks) * c!2 
+ etzs(jp,kp) * cl3 

else 

ell = delxi * delet 

clO = delxi - ell 

cOl = delet - ell 

cOO * 1 - delxi -delet + ell 

jp = jsp(js) 

kp = ksp(ks) 


xt 

= 

xyz ( js, l,Jcs, 1) 

★ 

cOO 

+ 

xyz ( jp. 

l,ks, l) 

★ 

clO 


+ 

xyz ( js, l,kp, 1) 

* 

cOl 

+ 

xyz ( jp. 

l»kp, 1) 

* 

ell 

yt 


xyz ( js, l,ks,2) 

★ 

cOO 

+ 

xyz ( jp. 

l,ks,2) 

* 

clO 


+ 

xyz ( js, l,kp,2) 

* 

cOl 

+ 

xyz ( jp. 

l,kp,2) 

* 

ell 

zt 

= 

xyz ( js, l,ks, 3) 

* 

cOO 

+ 

xyz ( jp. 

1, ks, 3) 

* 

clO 


+ 

xyz { js, l,kp, 3) 

* 

cOl 

+ 

xyz ( jp. 

1, kp, 3) 

* 

ell 


xix = xixs(js,ks) * cOO + xixs(jp,ks) * clO 

+ xixs(js,kp) * cOl + xixs(jp,kp) * ell 

xiy = xiys(js,ks) * cOO + xiys(jp,ks) * clO 

+ xiys(js,kp) * cOl + xiys(jp,kp) * ell 

xiz = xizs(js,ks) * cOO + xizs(jp,ks) * clO 

+ xizs(js,kp) * cOl + xizs(jp,kp) * ell 

etx = etxs(js,ks) * cOO + etxs(jp,ks) * clO 

+ etxs(js,kp) * cOl + etxs(jp,kp) * ell 

ety = etys(js,ks) * cOO + etys(jp,ks) * clO 

+ etys(js,kp) * cOl + etys(jp,kp) * ell 

etz = etzs(js,ks) * cOO + etzs(jp,ks) * clO 


1 



1 + etzs(js,kp) * cOl + etzs(jp,kp) * ell 

c endif 

ddelxi * xix * <xp - xt) 

1 + xiy * (yp - yt) 

2 + xiz * (zp - zt) 

ddelet = etx * (xp - xt) 

1 + ety * (yp - yt) 

2 + etz * (zp - zt) 

delxi * delxi + ddelxi 

delet * delet + ddelet 

c 

c Test for being far off. 
c 

if (max (abs (delxi-0 . 5) , abs (delet-0 . 5) ) .gt .3 . 5) goto 12 
c 

c Close enough? 
c 

error = ddelxi**2 + ddelet**2 
c error = sqrt (error) 

if (error . le . 1 .e-3 ) goto 12 

11 continue 

12 continue 
c 

c Converged to the right cell? 
c 

if (max (abs (delxi-0. 5) , abs (delet-0. 5) ) . le . 0 . 5+toler) goto 20 
c 

c Update cell number, 
c 

jso = js 
kso = ks 

if (delxi . It . 0 . ) then 

if (abs (delxi) .gt. toler) then 
js = jsr(js) 
endif 
endif 

if (delxi . gt . 1 . ) then 

if(delxi-l. .gt. toler) then 
js = jsp(js) 
endif 
endif 

if (delet . It . 0 . ) then 

if (abs (delet ) .gt. toler) then 
ks * ksr(ks) 
endif 
endif 

if (delet . gt . 1 . ) then 

if(delet-l. .gt. toler) then 
ks = ksp(ks) 
endif 
endif 
c 

c If (js r ks) didn't change, we hit the grid boundary. Limit delxi, delet 
c to the range [0,1]. 
c 

if ( js .eq. jso . and.ks .eq.kso) then 
delxi * min (max (delxi, 0 .), 1 . ) 
delet = min (max (delet , 0 . ) , 1 . ) 
istat = 1 

write (6, *) ' delxi, delet' , delxi, delet 
goto 20 
endif 


10 


continue 



write (6, *) 9 delxi, delet' ,delxi,delet 
c 

c Convergence failed. (Hope this never happens, but note it.) 
c 

ifail — ifail + 1 
if (ifail.gt.3) then 
istat * 2 
goto 20 
endif 

js = jstart - 1 
if(js .le. 0) js * 1 
ks * kstart - 1 
if (ks .le. 0) ks = 1 
goto 5 

20 continue 

c write{6,*) 9 js,ks' , js,ks 

c write (6,*) 'delxi delet' , delxi, delet 

return 
end 


subroutine heat (qw, qwbtu, wltemp, cf ) 
c 

c This program reads in streamline paths and metric coefficients, uses 
c Zoby' s approximate heating methodology and used Gupta's curve fit to 
c compute equilibrium air thermodynamics and transport properties to 
c evaluate convective heating along the streamline, 
c 

common /connty/ il,jl,kl 

common /stagpt/ xstag,ystag, zstag 

common /psidps/ psi<8), dpdxi(8), dpdeta<8), dpdzta{8) 
common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim,kdim, 

& iblank (151,40,110), blank, isubs (2) , 

& jsubs (2) , ksubs (2) 

common /qxyz/ u(151, 40,110) ,v(151,40,110) , 

& w (151, 40, 110) , p (151, 40, 110) 

common /save/ igsav, issav, isav, jsav, ksav, asav,bsav,gsav, 

& f irst, back, second, if atal 

common /freest/ pinf , t inf , rhoinf , vinf , vninf , cpinf , igas 
common /turbml/ iturb 
common /radeq/ itwall, wallt, irad, eps 
common /tempwl/ twall(1000) 
c 

c boundary and momentum thickness 

c 

common /bltrn/ deltal,thetl 
c 

dimension xs (1000) ,ys(1000) ,zs(1000) ,xn(1000),yn(1000) , 

& zn (1000) , vels (1000), ps (1000) ,hmet (1000) , 

& isaves (1000) , ksaves (1000) 

c 

logical blank, back, first , second 
c 

c conversion factors 
c 

data cmin, fps, atm, den /2.54, 30.48, 2116.2, .01601846/ 

data rtok /1.8/ 
c 

c Stef an-Bolt zman constant (w/cm**2/k) 
c 

data sigma /5.67e-12/ 
c 

rewind 15 
is = 1 



80 continue 

read(15,end~90) xs(is) ,ys(is) ,zs(is) , xn(is) ,yn(is) ,zn(is) , 

& vels (is) , ps (is) , hmet (is) , isaves (is) , ksaves (is) 

hmet ( is ) ~ hmet ( is ) *cmin 
is * is + 1 
go to 80 
90 continue 
it * 1 

92 continue 

if (itwall .eq. 1) then 

read (18, *, end=95) twall (it) 
twall(it) * twall (it) / rtok 
it « it + 1 
go to 92 
endif 

95 continue 

is = is -1 

if (itwall .eq. 0) then 
do 97 iu = l,is 
twall (iu) = wallt/rtok 
97 continue 

endif 
iwit = 0 
99 continue 

isubs(l) * 1 
isubs(2) * idim 
jsubs(l) * 1 
jsubs(2) * jdim 
ksubs(l) « 1 
ksubs(2) * kdim 
blank = .false, 
c 

igsav - 1 
issav = 1 
asav = 0 . 
bsav = 0 . 
gsav * 0 . 
s = 0 . 
sumO = 0 . 
fl = 0. 
ite = 0 
iweq 38 0 

do 200 iter = l,is 
rewind 6 

write (6,*) 9 iter = 9 , iter 
icon *= 0 
c 

c initial edge velocity 
c 

vele = vels (iter) * vninf 
velet * vele 

c vele = vels (iter) * 30.48 

c 

c edge pressure 
c 

pe = ps(iter) * pinf 
c 

c stagnation enthalpy (cal/g) 
c 

hinf = cpinf*tinf + 0 . 5*vinf **2*2 . 3901e-8 
write (6,*) ' hinf (computed) = ',hinf 
c 

c initial search cell 
c 

isav = isaves (iter) 
jsav = 1 



ksav « ksaves(iter) 

100 continue 
c 

c edge enthalpy (cal/g) 
c 

c he = hinf - 0 . 5*vele**2*2 . 3901e-8 

he = hinf - 0 . 5*velet**2*2 . 3901e-8 
c 

c get edge property 
c 

he * 1.0e-3*he 
write (6,*) ' he * ',he 
iflag * 1 

if(igas .eq. 0) then 
iflag * 0 
te * he /.24e-3 
else if (he .It .12) then 
te * he / .24e-3 
else 

call enthalpy (pe, he, te, if lag) 
if(ifatal .ne. 0) return 
endif 

call eqprop (pe, te, cpe, rhoe, xmue, xke, pre, gammae, airme, iflag) 
if(ifatal .ne. 0) return 
c 

c get wall properties (tw is input value) 
c 

iflag * 2 

if(igas .eq. 0) iflag * 0 
tw * twall(iter) 
call enthalpy (pe, hw, tw, if lag) 
if(ifatal .ne. 0) return 

call eqprop (pe, tw, cpw, rhow, xmuw, xkw, prw, gammaw, airmw, iflag) 
if(ifatal ,ne. 0) return 
c 

c compute Eckert's reference enthalpy 
c 

c hstar - 0.5*(he+hw) + . ll*sqrt (prw) *vele**2 * 2.3901e-ll 

hstar = 0.5*(he+hw) + . ll*sqrt (prw) *velet**2 * 2.3901e-ll 
c 

c get reference properties 
c 

iflag = 1 

if(igas .eq. 0) then 
tstar = hstar / .24e-3 
iflag * 0 
else 

call enthalpy (pe, hstar, tstar, iflag) 
if(ifatal .ne. 0) return 
endif 

call eqprop (pe, tstar, cpstar, rhostr, xmustr, xkstar,prstar,gamstr , 
& armstr, if lag) 

if (ifatal .ne. 0) return 
c 

c compute momentum boundary layer thickness 
c 

if (iter .eq. 1) then 

ds = sqrt ( (xstag-xs (iter) ) **2 + (ystag-ys (iter) ) **2 
& + (zstag-zs (iter) ) **2 ) 

ds = ds*cmin 

s = s + ds 

c sum = 0.5*f2*ds 

c sum = 0.25*f2*ds 

else 

ds = sqrt ( (xs (iter) -xs (iter-1) ) **2 + (ys (iter) -ys (iter-1) ) **2 

& + (zs (iter) -zs (iter-1) ) **2) 
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ds - ds*cmin 
s ~ s + ds 

c sum * sumO + 0 . 5* (f 1+f 2) *ds 

endif 
c 

c check if turbulence heating is required 
c 

rein * 1 . 
c3 * 2. 
c4 - 0.5 
c2 * .440896 
iterb = 0 
thetlo - 0 . 

110 continue 

c f2 = rhostr*xmustr**rem*vele*hmet (iter) **c3 

f2 * rhostr*xmustr**rem*velet*hmet (iter) **c3 
if (iter .eq. 1) then 
c sum * 0.25*f2*ds 

sum = 0.5*f2*ds 
else 

sum * 0.5 * (fl + f2) * ds + sumO 
endif 

thetl = (c2*sum) **c4/ (rhoe*vele*hmet (iter) ) * .93 
thetl = (c2*sum) **c4/ (rhoe*vele*hmet (iter) ) 
thetl = (c2*sum) **c4/ (rhoe*velet*hmet (iter) ) 
write (6,*) ' thetlo, thetl « ', thetlo, thetl 
reth = rhoe * vele * thetl / xmue 
reth — rhoe * velet * thetl / xmue 
if(iturb .ne. 0) then 

if (abs (thetl-thetlo) .gt. .01 .or. iterb .eq. 0) then 
rncap - 12 . 67-6 . 5*alogl0 (reth) +1 .21* (aloglO (reth) ) **2 
rncap *» anint (rncap) 
rem * 2. / (rncap + 1.) 
c5 * 2.2433 + 0.93 * rncap 
c3 ■ 1 . + rem 
c4 *= 1 . / c3 

cl = { (1 . /c5) ** (2 . *rncap/ (rncap+1 . ) ) ) * (rncap/ ( (rncap+1 . ) * 

& (rncap+2 . ) ) ) **rem 

c2 = (l.+rem)*cl 
iterb = iterb + 1 
write (6,*) 9 iterb' , iterb 

write (6 , *) 9 rncap, rem, cl, c2, c3, c4, c5' , rncap, rem, cl, c2, c3, c4,c5 
if (iterb .ge. 20) then 

write (6,*) 9 Too many iterations !!!' 
ifatal = 1 

if(iwit .eq. 0) then 
qw = 0 . 
qwbtu = 0 . 
wltemp = 0 . 
cf - 0. 
endif 
return 
endif 

thetlo = thetl 
go to 110 
endif 
endif 

f2 = rhostr*xmustr*vele*hmet (iter) **2 
thetl - . 664*sqrt (sum) / (rhoe*vele*hmet (iter) ) * . 93 
thetl = . 664*sqrt (sum) / (rhoe*vele*hmet (iter) ) 
write (6,*) ' Momentum thickness » thetl 

confute boundary layer thickness 

if(iturb .eq. 0) then 
deltal = 5.55*thetl 



else 

c haw * he+O . 5*sqrt (prw) *vele**2*2 . 3901e-ll 

haw — he+O . 5*sqrt (prw) *velet ** 2 * 2 . 3901e-ll 
thedd * rncap + 1. + ( ( (rncap+2 . ) *hw/ (rncap*haw) +1 . ) * 
c & (1 • +1 . 29*prw** . 333* . 5*vele**2*2 . 3901e-ll/he) ) 

& (1 . +1 . 29*prw** . 333* . 5*velet**2*2 . 3901e-ll/he) ) 

deltal * thetl * thedd 
endif 

c deltal = 5 . 55*thetl*l . 18 

write (6,*) ' Boundary layer thickness = ', deltal 
c deltal * 0 . 

if (icon .eq.l) go to 150 
c 

xble - xs(iter) + xn (iter) *deltal/cmin 
yble = ys(iter) + yn (iter) *deltal/ciain 
zble = zs(iter) + zn (iter) *deltal/cmin 
c 

c convert back to inch for interpolation 
c 

c xble = xble/cmin 

c yble = yble/cmin 

c zble = zble/cmin 

c 

c find cell for interpolation 
c 

isrch = 0 
igrid = 1 
nsubs = 1 

call close3 (isrch, igrid, nsubs, i, j , k, a, b, g, xble, yble, zble, isub, 

& istat) 

if(istat ,ne. 0) then 

write (6,*) ' close3 search failed' 
c stop 

i = isaves(iter) 

j - jj 

c j * jdim - 1 

if ( j .eq. 0) j - 1 
k = ksaves (iter) 

if (ksaves (iter) .eq. kdim) k - k - 1 
write(6,*) ' i,j,k = f ,i,j,k 
else 

write ( 6, *) ' i,j,k = ',i,j,k 

jj = j 

endif 

c 

c interpolation of edge properties 
c 

il = i 
jl - j 

kl = k 

call find (xble, yble, zble) 
if(ifatal .ne. 0) return 
c velb = vele 

velb = velet 
ipl = il + 1 

jpl = jl + 1 

kpl = kl + 1 

ue = u (il, jl, kl) *psi (1) +u (ipl, jl,kl) *psi (2) +u (ipl, jl, kpl ) *psi (3) 

& + u (il, jl,kpl) *psi (4) +u (il, jpl, kl) *psi (5) +u (ipl, jpl, kl) *psi (6) 

& + u (ipl, jpl, kpl) *psi (7) +u (il, jpl, kpl) *psi (8) 

c 

ve = v (il, jl, kl) *psi (1) +v(ipl, jl,kl) *psi(2) +v(ipl, jl,kpl) *psi (3) 

& + v(il, jl, kpl ) *psi (4) +v (il, jpl,kl) *psi (5) +v(ipl, jpl,kl) *psi (6) 

& + v (ipl, jpl, kpl) *psi (7) +v (il, jpl, kpl) *psi (8) 

c 

we = w (il, jl, kl ) *psi (1) +w (ipl, jl, kl) *psi (2 ) +w (ipl, jl,kpl) *psi(3) 
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& + w(il f jl, kpl) *psi (4) +w(il, jpl,kl) *psi (5)+w{ipl r jpl/kl) *psi (6) 

& + w(ipl, jpl/kpl) *psi (7) +w(il, jpl, kpl) *psi (8) 

c 

pe = p (il/ jl, kl) *psi <1) +p (ipl/ jl/ kl) *psi (2) +p (ipl, jl, kpl) *psi (3) 

4 + p (il, jl, kpl) *psi (4) +p (il, jpl, kl) *psi (5) +p (ipl/ jpl f kl) *psi (6) 

+ p (ipl, jpl, kpl) *psi (7) +p (il, jpl/kpl) *psi (8) 
c 

pe = pe*pinf 
c 

vele = sqrt(ue*ue + ve*ve + we*we) 
c 

c use tangential component 
c 

velen = ue*xn(iter) + ve*yn(iter) + we*zn(iter) 
velet = sqrt (vele*vele - velen*velen) 
velet = velet *vninf 
vele * vele*vninf 
c vele *= vele*30.48 

c 

write (6/*) 9 vele, velb = ', vele, velb 
c tconv = (vele - velb) /velb 

tconv = (velet - velb) /velb 
if (abs (tconv) .It. .001) then 
icon * 1 
else 

ite = ite + 1 
if(ite .gt. 25) then 
c write (6/*) 9 edge failed' 

icon = 1 
endif 
endif 
go to 100 
c 

c compute laminar heating 
c 

150 continue 
c 

c momentum thickness Reynolds number 
c 

c reth = rhoe*vele*thetl/xmue 

write (6,*) 9 reth = ' / reth 
c 

compute internal energy per unit mass (m**2/s**2) 

ein - he*4 . 184el0 - pe*l . 01336e6/rhoe 
eins = ein * l.e-4 
rhoes = rhoe * l.e3 
c 

c use tannelhill's curve fit to evaluate entropy 
c 

call tgas2 (eins / rhoes/ entrpy) 
c 

c adiabatic wall enthalpy (kcal/g) 
c 

c haw = he +0 . 5*sqrt (prw) *vele**2*2 . 3901e-ll 

haw = he +0 . 5*sqrt (prw) *velet**2*2 . 3901e-ll 
if (iturb .eq. 0) then 

cf = .44 / reth* (rhostr/rhoe) * (xmustr/xmue) 
if (xs (iter) .le. 128.8) then 

cf = 1 . 534* (rhoe*xmue) ** . 43* (xmuw*rhow) ** . 07 
endif 

c cf = .44 / reth 

c qw = 0 . 22 * rhost r *xmust r * rhoe*vele/ (rhoe*xmue) * 

qw = 0 . 22*rhostr*xmustr*rhoe*velet / (rhoe*xmue) * 

& (haw - hw) / (reth*prw** . 6) 

else 



cf = 2 . *cl/ (reth**rem) * (rhostr/rhoe) * (xmustr/xmue) **rem 
c cf = 2 . *cl/ (reth**rem) 

c qw = cl* (rhostr/rhoe) * (xmustr/xmue) **rem*rhoe*vele* 

qw = cl* (rhostr/rhoe) * (xmustr/xmue) **rem*rhoe*velet* 

& (haw-hw) / (prw** . 4*reth**rem) 

endif 
c 

c convert kcal/cm**2/s to watt/cm**2 
c 

qw = qw*4.18392e3 
c 

c confute radiation equilibrium wall temperature 
c 

if (qw .le. 0.) then 
tweq = 300 . 
qw = 0. 
else 

tweq *= <qw/ (eps*sigma) ) ** .25 
endif 

if (abs (tweq-twall (iter) ) ,gt. 5.) iweq = 1 
twall(iter) = tweq 
write (25,*) tweq 
c 

c convert to btu/ft**2/s 
c 

qwbtu = qw* . 88 
s = s + ds 
sumO = sum 
fl = f 2 

velets *= velet 
ite = 0 

wltemp = twall(iter) 
fss = xmue/xmustr 
rethss = fss*reth 
rethl = aloglO (rethss) 

c cf c = (te/tstar) / (17 . 08*rethl*rethl + 25.11*rethl + 6.012) 

write (24,9001) xs (iter) , s, qw, qwbtu, cf 
write (34,9002) xs (iter) ,ys (iter) , zs (iter) ,entrpy 
200 continue 

if(irad .ne. 0) then 
if (iweq .ne. 0) then 
iweq = 0 
iwit ~ iwit + 1 
if(iwit .gt. 7) then 

write (6,*) 9 Wall equilibrium temperature iteration failed' 
go to 9999 
endif 
rewind 25 
rewind 24 
rewind 34 
go to 99 
endif 
endif 
c 

9001 format (5 (lx, el 4 . 7) ) 

9002 f ormat ( 4 (lx, el4 . 7) ) 

9999 continue 

return 

end 

c 

SUBROUTINE TGAS2(E,R,S) 

C 

C INPUTS FOR SUBROUTINE: 

C E=INTERNAL ENERGY IN (M/SEC) **2 

C R=DENSITY IN KG/M* *3 

C 



c 

c 

c 

c 


10 


20 


30 

40 


50 


60 

70 


80 


90 


100 


OUTPUT: 

S “ENTROPY IN (M/SEC) **2/K 


DATA EO, R0, GASCON/7 8408 .4E00, 1.2 92E00, 287. 06E00/ 

RRATIO=R/R0 

ERATIO-E/EO 

Y-ALOG10 (RRATIO) 

Z»ALOG10 (ERATIO) 

IF (ABS (Y+4 . 5E00) .LT.2.5E-02) GO TO 10 

IF (ABS (Y+0 . 5E00) .LT.0.5E-02) GO TO 40 

IFLAG=-1 

GO TO 80 

IFLAG-0 

RSAVE=R 

YM=Y 

Y— 4 . 5E00+2 . 5E-02 

YHIGH=Y 

R= (10 . **Y) *R0 

JFLAG=-1 

GO TO 80 

SHIGH=S 

Y— 4.5E00-2.5E-02 
YLOW=Y 

R= (10 . **Y) *R0 

JFLAG“0 

GO TO 80 

SLOW=S 

GO TO 70 

IFLAG=1 

RSAVE=R 

YM-Y 

Y— O.5EOO+0 .5E-02 

YHIGH=Y 

R= (10 . **Y) *R0 

JFLAG—1 

GO TO 80 

SHIGH=S 

Y=-0.5E00-0 .5E-02 
YLOW-Y 

R=(10.**Y) *R0 
JFLAG=0 
GO TO 80 
SLOW=S 

S=SLOW+ ( SHIGH— SLOW) / (YHIGH-YLOW) * (YM-YLOW) 

R=RSAVE 

RETURN 

CONTINUE 

IF (Z .LE . 0 . 65E00) GO TO 110 
IF (Y.GT .-4 . 5E00) GO TO 90 
IF (Z.GT . 3 . 69E00) WRITE ( 6, 1000) R,E 
GASl - — 9 . 91081E-01-5 . 00277E00*Y 
GAS2* (5 . 46521E01+5 . 10144E00*Y) *Z 

GAS3= (1 . 76206E-02+2 . 12002E-02*Z+1 . 76358E-03*Y) *Y*Y 
GAS4= (-2 . 97001E01-1 . 84915E00*Y+5 . 87892E00*Z) *Z*Z 
GO TO 120 

IF (Y.GT .-0 . 50EOO ) GO TO 100 
IF (Z .GT . 3 . 4E00) WRITE (6 f 1000) R,E 
GAS1>=1 . 083 6E0 1-4 . 5552 4E00^Y 
GAS2=(2.96473E01+3.90851E00*Y) *Z 

GAS3= (-2 . 05732E-03+3 . 65982E-02*Z+5 . 23821E-03*Y) *Y*Y 
GAS4= (-1. 6700 1E0 1-1. 44623E00*Y+3 .98307EOO*Z) *Z*Z 
GO TO 120 

IF (Z .GT . 3 . 0E00) WRITE (6, 1000) R,E 
GAS1=2 . 01858E01-3 .13458E00*Y 
GAS2= (1. 0361 9E01+1. 877 67E00*Y) *Z 



GAS 3= {-1 . 72922E-01+1 . 12174E-01*Z+1 . 28 62 6E-02*Y) *Y*Y 
GAS 4= (-5 . 43557E00-8 . 71048E-01*Y+2 . 01789EOO*Z) *Z*Z 
GO TO 120 

110 DELTZ=Z-0 . 4E00 

DELTS= (2 . 5EOO*DELTZ-Y) *GASCON*2 . 302585E00 
S=677 9 . 2004E00+DELTS 
GO TO 130 

120 SNON=GAS 1 +GAS 2 +GAS 3 4-GAS 4 

S=GASCON* SNON 

130 IF (IFLAG) 160, 140, 150 

140 IF(JFLAG) 20,30,160 

150 IF ( JFLAG) 50,60,160 

160 CONTINUE 

1000 FORMAT </20X, 48HWARNING] OUTSIDE VALIDITY RANGE OF CURVE FIT 
*, /,20X,5HRHO =, 1PE15 . 8, 5X, 3HE =,1PE15.8,/) 

RETURN 

END 


c 

subroutine enthalpy (pinp, hinp, t , if lag) 
c 

c This subroutine takes input pressure and enthalpy, uses Gupta's 
c equilibrium air thermodynamic and transport properties curve fit 
c to compute temperature, 
c if lag = 0, perfect gas 

c iflag = 1, pressure and enthalpy are input 
c iflag - 2, pressure and temperature are input 
c 

dimension ptb(7), ah (6, 7), bh(6,7), ch(6,7), dh(6,7), eh (6, 7), 
& temp(6,7), tmax(6,7), hmax<6,7), ntmax(7) 

c 

common /save/ igsav, issav, isav, jsav, ksav, asav,bsav,gsav, 

& first, back, second, if atal 

common /coef/ coefl, coef2, coef3, coef4, coefS 
common /pass/ ipass 


c 


data 

ptb / 

rH 

o 

o 

o 

.001, . 

01, .1, 

1., 10., 100./ 

data 

ntmax 

16 , 6, 

5, 4 , 

4, 4, 3/ 


data 

ipmax 

m 




data 

tmax 

/2250 . , 

4250., 

6750., 

10750., 17750., 25000., 

& 


2250., 

4250., 

6750. , 

11750., 18750., 28000., 

& 


2750., 

5250., 

9750., 

17750., 30000., 30001., 

& 


3250., 

6250., 

15250., 

30000., 30001., 30001. 

& 


3750., 

8250., 

17750., 

30000., 30001., 30001. 

& 


4250., 

9250., 

18750., 

30000., 30001., 30001. 

& 


6250., 

12750. 

, 30000. 

, 30001., 30001., 30001 


data ah /.128180el, .125380e2, .426138e2, .885088el, .151569e2, 

& . 101759e2, 

& .902850, . 237222e2, .880011e2, -.333238e2, .196866e2, 

& . 446849e2, 

& .653358, . 431122el, -.126229el, .209845e2, .268647e2, 0 

& .363885, - . 8 65884el, -.164319e2, -.207249e2, 0., 0., 

& .209284, - . 171560e2, -.134978e2, -.564265el, 0., 0., 

& .124937, - . 120314e2 , -.913636el, .639208el, 0., 0., 

& - . 755123e-2, -.117469el, -.245329el, 0., 0., 0./ 


data 

& 

& 

& 

& 

& 0 ., 

& . 329839el, 

& . 187458el, 

& . 109286el, 

& . 164258e-l 


. 720107e2, .123001e3, 

. 118014e3, . 213329e3, 

.267 604e2, .113432e2, 

- . 208034e2, -.285858, 

-. 416138e2, .801118el, 

- . 229170e2, .113996e2, 

- . 592622el, .371340el, 


- . 207380e2, -.713138el 

- . 316397e2, -.201771e2 

- . 181381e2, - . 104256e3 

. 633182e2, 0., 0., 

. 262889e2, 0., 0., 

- . 149544e2 , 0., 0., 
0 ., 0 ., 0 ./ 


bh / . 121182e2, 

- . 161 956e2, 
. 839944el, 

- . 141086e3, 
. 59688 6el, 



c 


c 


c 


c 


data ch 

& 

& 

& 

6 

& 

& 

& 

& 

& 


/ . 424907e2, .148949e3, .121801e3, -.13460402, -.172524, 
-.336892el, 

. 289458e2 f .214780e3, .181623e3, -.401000el, .635249el, 
. 159412e3, 

.201689e2, .541203e2, .109117e2, -.399635, .145439e3, 

0 ., 

. 110641e2, - . 132700e2, .447878el, -.678713e2, 0., 0., 

. 622153el, -.332532e2, .192371el, -.396119e2, 0., 0., 

. 355163el, -.129249e2, -.259796el, .882252el, 0., 0., 
.366590, - . 214181el, -.288683, 0., 0., 0./ 


data dh 

/ . 665524e2, 

. 133853e3, 

. 509305e2, 

. 166408el, .643645, 

& 

& 

. 161274e2, 
. 448640e2, 

. 171168e3, 

. 661367e2, 

. 379639el, -.174347, 

& 

& 

- . 738595e2, 
. 309518e2, 

. 462077e2, 

. 400303el, 

. 387388el, -.846045e2 

& 

& 

0 . , 

. 173605e2, 

.242899el, 

. 196275el, 

. 312942e2, 0., 0., 

& 

. 101561e2 , 

- . 747816el 

, .930272, 

. 251297e2, 0., 0., 

& 

. 617946el, 

.262066, . 

114665el, . 

258596el, 0., 0., 

& 

. 210603el, 

.251111el, 

.421200, 0 

., 0., 0./ 

data eh 

/ . 385195e2, 

. 451550e2 , 

. 995964el, 

. 356570el, .356353el, 

& 

& 

- . 201068el, 
. 256452e2 , 

. 513939e2, 

. 110476e2, 

. 3254 69el, .354258el, 

& 

& 

r 

. 155141e2, 
. 174843e2, 

. 152182e2, 

.284253el, 

.283981el, .212051e2, 

« 

& 

U . , 

. 999025el, 

. 417259el, 

.256061el, 

- . 158288el, 0., 0., 

& 

. 603650el, 

. 178858el, 

. 244209el, 

-.207198el, 0., 0., 

& 

. 386028el, 

.235363el, 

.236890el, 

. 107086el, 0., 0., 

& 

. 195195el, 

.212013el, 

. 239842el, 

o 

Si 

• 

O 

o 

pinps * 
if (if lag 

0. 

.eq. 1) write (6,*) r pinp, hinp = 

' , pinp, hinp 


if (if lag .eq. 0) then 
hinp * . 24e-3*t 
return 
endif 


c 

c find pressure range for interpolation 
c 

if (pinp .It. l.e-4) then 
pinps = pinp 
pinp = l.e-4 
endif 

if (pinp .It. l.e-4) then 

if (if lag .eq. 2 .and. t ,le. 500.) then 
hinp = ,24e-3*t 
else 

write (6,*) 9 Pressure is below the lower bound, p — ',pinp 
ifatal = 1 
return 
endif 

else if (pinp .gt. 100.) then 

write (6,*) 9 Pressure is above the upper bound, p = ',pinp 
ifatal = 1 
return 
endif 


c 

c find the index which brace the input pressure 
c 

do 10 i = 1,7 

if (pinp ,le. ptb(i)) then 

if (abs (pinp-ptb <i) ) .le. l.e-4) then 



intp = 0 
else 

intp * 1 
endif 

ip = i 

ipml ■* ip - 1 
go to 20 
endif 

10 continue 

write (6,*) 9 Pressure is out of range' 

20 continue 

c write (6,*) ' ip, ipml ', ip, ipml 

if(iflag .eq. 2) then 
xi = alog(t* . 0001) 
if (intp .eq. 0) go to 50 
do 30 i «= 1 , ntmax ( ipml ) 
if (t .gt. tmax (i, ipml) ) go to 30 
it = i 
go to 40 
30 continue 

write (6,*) 9 Temperature ',t,' is outside the range of', 

&' available data' 
ifatal = 1 
return 

40 continue 

if (t .It. 500.) then 
hinp = .24e-3*t 
if (pinps ,ne. 0.) pinp = pinps 
return 
endif 

hi = ah (it, ipml) *xi**4 + bh (it, ipml) *xi**3 + ch (it, ipml)* 

& xi**2 + dh (it, ipml) *xi + eh (it, ipml) 

c 

50 continue 

if (t .It. 500.) then 
hinp — . 24e-3*t 
if (pinps ,ne. 0.) pinp = pinps 
return 
endif 

do 60 i = 1, ntmax (ip) 
if (t .gt. tmax (i, ip)) go to 60 
it = i 
go to 70 
60 continue 

write(6,*) ' Temperature ',t,' is outside the range of', 

&' available data' 
ifatal = 1 
return 

70 continue 

h2 = ah (it , ip) *xi**4 + bh (it , ip) *xi**3 + ch (it, ip) *xi**2 
& + dh(it,ip)*xi + eh (it, ip) 

c 

if (intp .eq. 0) then 
alogh = h2 
else 

alogh = (h2 - hi) * (alog(pinp) - alog (ptb (ipml ) ) ) 

& / (alog (ptb (ip) ) - alog (ptb (ipml) ) ) + hi 

endif 
c 

hinp = exp (alogh) 
go to 200 
endif 
c 

c compute maximum enthalpy for each temperature range for each pressure 
c 


C'Z 


if(ipass .eq. 0) then 



ipass = 1 

do 80 j * l,ipmax 

do 80 i * l,ntmax(j) 

xi * alog(tmax(i, j) *.0001) 

hmax(i,j) - exp (ah (i f j) *xi**4 + bh (i, j) *xi**3 + 

& ch(i, j) *xi**2 + dh (i, j) *xi + eh(i, j)) 

c write(6,*) 9 hmax 9 , j, i, hmax (i, j) 

80 continue 

endif 

if(intp .eq. 0) go to 100 
c 

c find correct temperature range and curve to solve for temperature 
c 

ntip = ntmax(ipml) 
do 90 i * l,ntip 
if(hinp .It. hmax(i, ipml) ) then 
coefl * ah (i, ipml) 
coef2 * bh(i f ipml) 
coef3 = ch(i,ipml) 
coef4 * dh(i,ipml) 
coef5 = eh (i, ipml) 
tmin = 500 . 
hmin * .24e-3*tmin 
if (i .gt. 1) then 

tmin * tmax (i-1, ipml) 
hmin = hmax (i-1 , ipml) 
endif 

c ti * (tmax (i, ipml) +tmin) * . 5 

ti=tmin+ (a log (hinp) -a log (hmin) ) * (tmax (i, ipml) -tmin) 

& / {alog (hmax (i, ipml) ) -alog (hmin) ) 

call newton (ti, hinp, tl) 
go to 100 

else if (hinp .eq. hmax ( i , ipml ) ) then 
tl = tmax (i, ipml) 
go to 100 
endif 

90 continue 
c 

write (6,*) ' Enthalpy is out of range', hinp 

ifatal * 1 

return 

100 continue 
c 

ntip = ntmax(ip) 
do 110 i = l,ntip 
if (hinp .It. hmax (i, ip)) then 
coefl = ah(i,ip) 
coef2 = bh(i,ip) 
coef3 = ch(i,ip) 
coef4 = dh(i,ip) 
coefS = eh(i,ip) 
tmin = 500. 
hmin * .24e-3*tmin 
if(i .gt. 1) then 
tmin = tmax(i-l,ip) 
hmin = hmax(i-l,ip) 
endif 

ti *= tmin + (alog (hinp) -alog (hmin) )* (tmax (i, ip) -tmin) 

& / (alog (hmax (i, ip) ) -alog (hmin) ) 

call newton (ti, hinp, t2) 
if (ifatal .ne. 0) return 
go to 120 

else if (hinp .eq. hmax (i, ip)) then 
t2 = tmax (i, ip) 
go to 120 
endif 


110 continue 
c 

write (6,*) ' Enthalpy is out of range' , hinp 

ifatal = 1 

return 

120 continue 

if(intp .eq. 0) then 
t = t2 
go to 200 
endif 
c 

c use logrithm interpolation to find temperature 
c 

alogt - (alog(t2) - alog(tl)) * (alog(pinp) - alog (ptb <ipml) ) ) 

& / (alog (ptb (ip) ) - alog (ptb (ipml) ) ) + alog(tl) 

c 

t = exp (alogt) 
c 

200 continue 
c 

if (pinps .ne. 0.) pinp = pinps 
return 
end 
c 

subroutine newton (ti, hinp, t) 

common /save/ igsav, issav, isav, jsav, ksav, asav,bsav, gsav, 

& first, back, second, if atal 

c 

c use newton's method to solve for polynomial 
c 

common /coef/ a, b, c, d, e 
c 

xi = alog(ti* . 0001) 
iter = 0 
10 continue 

resl = exp(a*xi**4 + b*xi**3 + c*xi**2 + d*xi + e) 
res = resl - hinp 

if(abs(res) .It. 1.0e-4) go to 100 
c 

c compute derivative 
c 

fp = resl * (4 . *a*xi**3 + 3.*b*xi**2 + 2.*c*xi + d) 

dxi = -res/fp 

xi = xi + dxi 

iter » iter + 1 

if (iter .gt. 50) then 

write (6,*) 9 Newton method failed to converge' 
ifatal = 1 
return 
endif 
go to 10 
100 continue 

t = 10000. * exp(xi) 

return 

end 

subroutine eqprop (pinp, tinp, cp, rho, xmu, xk,pr, gamma, airm, if lag) 
c 

c This routine uses Gupta's curve fit equation to obtain thermodynamics 
c and transport properties using logarithm interpolation, 
c Inputs are pressure (pinp) and temperature (tinp) 
c 

common /save/ igsav, issav, isav, jsav, ksav, asav, bsav, gsav, 

& first, back, second, ifatal 

c 

dimension acp(9,7), bcp(9,7), ccp(9,7), dcp(9,7), ecp(9,7), 

& az(5,7), bz(5,7), cz{5,7), dz(5,7), ez(5,7), fmu(4,7). 



*> fh X> 


amu(4,7), bmu(4,7), emu (4, 7), dmu{4,7), emu(4,7), 
ak { 7 , 7 ) , bk (7,7) , ck<7,7), dk<7,7), ek<7,7), 
apr (8, 7) ,bpr (8,7) , cpr<8,7) ,dpr (8,7) ,epr<8,7) , fpr (8, 7) , 
ptb(7), tcpmax{9, 7) , tzmax(5,7), tmumax ( 4 , 7 ) , 

& tkmax (7, 7), tprmax(8,7), ncpmax(7), nzmax(7), 

& nmumax ( 7 ) , nkmax ( 7 ) , nprroax { 7 ) 

c 

c runiv is in cal/g-mole-k 

c 

data airmO, runiv /28.96, 1.987/ 

data ptb /l.e-4, l.e-3, l.e-2, l.e-1, l. f 10./ 100./ 
data ncpmax /9, 8 , 1, 8/ 1, 7, 6/ 
data nzxnax /5, 5/ 5 , 5, 5, 4, 3/ 
data nmumax /4 r 4/ 4, 4/ 3 f 3, 2/ 
data nlanax /l, 1 , l t € r 6 f 5/ 4/ 
data nprmax /8, 8/ S f 1, 1, 5, 5/ 
data ipmax 111 
c 

data tcpmax / 

& 12 50. /1750. ,2750. ,4750. , 6250. , 9750. ,142 50. /I 97 50. ,25000., 

& 1250., 2250., 3750., 5250., 7250., 10750., 17250. ,28000., 28001., 

& 1750., 2750., 4750., 6750., 12750., 19750., 30000., 30001., 30001. , 

& 1750., 2750., 4250., 6750., 9750., 15750., 2 1500., 30000., 30001., 

& 1750., 3250., 4750., 7750., 11750., 20500., 30000., 30001., 30 001., 

& 1750., 3250., 5750., 9250., 13750., 22500., 30000., 30001., 30001., 

& 1750., 3750., 6750., 10750., 17750., 30000., 30001., 30001., 30001./ 
c 

data tzmax / 

& 2750., 5750., 8750., 17750., 25000., 

& 3250., 6750., 9750., 19750., 28000., 

& 3250., 7250., 11750., 21500., 30000., 

& 3750., 8250., 13750., 23500., 30000., 

& 5750., 9250., 15750., 23500., 30000., 

& 5750., 9750., 17250., 30000., 30001., 

& 8750., 17750., 30000., 30001., 30001./ 
c 

data tmumax / 

& 7750., 10750., 16750., 25000., 

& 8250., 12250., 18750., 28000., 

& 8750., 14250., 19750., 30000., 

& 9750., 16750., 24500., 30000., 

& 11250., 19750., 30000., 30001., 

& 12750., 21500., 30000., 30001., 

& 15250., 30000., 30001., 30001./ 
c 

data tkmax / 

& 1750., 2750., 4750., 6250., 10250., 17750., 25000., 

& 1750., 2750., 4750., 6250., 11250., 18250., 28000., 

& 2250., 3250., 5750., 7750., 12750., 18750., 30000., 

& 2250., 4250., 6750., 9250., 16750., 30000., 30001., 

& 2250., 4250., 7750., 10750., 19250., 30000., 30001., 

& 3250., 5250., 8750., 13750., 30000., 30001., 30001., 

& 3750., 6250., 10750., 30000., 30001., 30001 ., 30001 . / 
c 

data tprmax / 

& 2250., 3750., 5750., 8250., 10750., 14750., 18250., 25000., 

& 2250., 4750., 7250., 10250., 12750., 17250., 20500., 28000., 

& 2750., 5250., 8250., 11750., 14250., 18250., 23500., 30000., 

& 2750., 5250., 7750., 13750., 18250., 25500., 30000., 30001., 

& 2750., 4750., 7750., 13250., 17750., 23500., 30000., 30001., 

& 2750., 5750., 10750., 20500., 30000., 30001., 30001., 30001., 

& 2750., 6750., 12750., 20500., 30000., 30001., 30001., 30001./ 
c 

data acp / 

& .349023, . 152264e2, -.159675e2, -.108293e3, -.116264e4, 

& - .238707e2, -.209557e2, .762671e3, -.789820e3, 


ft* 


4 .199532, .345376el, -.369572e2, -.146237e3, -.758521e3, 

& ~.330240e2, -.618098e2, .125063e3, 0., 

4 .669436, - . 453138e2, -.151035e3, .539167e3, .217707e2, 

4 -.122810e3, .162348e3, 0., 0., 

4 .291577, - . 662937el, .128388e3, -.296048e2, -.308894e3, 

4 .104767e3, -.188079e3, .232697e3, 0., 

4 .164992, - . 830572el, .848335e2, -.945467el, -.153176e3, 

4 . 975058e2, -.473648e2, 0., 0., 

& .111751, .252675, .450386e2, .231376e2, -.799940e2, 

& . 491689e2, -.253231e3, 0., 0., 

& . 986591e-l, .974261e-l, .210207e2, .143729e2, -.347606e2, 

4 . 450529e2, 0., 0., 0./ 

data bcp / 

& .344158el, .129277e3, -.136508e3, -.515276e3, -.266973e4, 
4 -.104336e3, .253228e2, -.167407e4, ,263864e4, 

4 . 192597el, .315624e2, -.128366e3, -.581296e3, -.139794e4, 
4 -.866157e2, .103127e3, -.298121e3, 0., 

4 . 644478el, -.292666e3, -.591051e3, .126894e4, -.450370e2, 

4 . 240030e3, -.497482e3, 0., 0., 

4 .278787el, -.382984e2, .596922e3, -.133243e3, -.267701e3, 
& - . 105447e3, .472158e3, -.869061e3, 0., 

& . 156336el, -.483112e2, .361629e3, -.640807e2, -.476111e2, 
& - . 158721e3, .818135e2, 0., 0., 

& . 105018el, . 341131el, .167261e3, -.104484el, .170114e2, 

& - . 116351e3, . 955890e3, 0., 0., 

.923581, . 14677 6el, .677318e2, -.128820e2, .320177e2, 

-. 143364e3, 0., 0., 0./ 

data ccp / 

& . 126715e2, .411057e3, -.411657e3, -.882748e3, -.221802e4, 
& -.890658e2, .212355e2, .130713e4, -.326378e4, 

& . 694347el, .107177e3, -.129698e3, -.848597e3, -.900003e3, 
& -.489572e2, -.262275e2, .210795e3, 0., 

& . 230631e2, -.699603e3, -.835692e3, .106221e4, -.192634e2, 
& -.138486e3, .525270e3, 0., 0., 

& . 992221el, -.779456e2, .101945e4, -.187832e3, -.478605e2, 
4 . 127166e2, -.407311e3, .117775e4, 0., 

4 . 552429el, -.101598e3, .561712e3, -.893740e2, .217674e2, 

4 .753693e2, .169726e2, 0., 0., 

4 .368846el, .131529e2, .224425e3, -.271807e2, .187072e2, 

4 .889977e2, -.132457e4, 0., 0., 

4 . 323392el, .575473el, .778089e2, -.173603e2, .148249el, 

4 . 161302e3, 0., 0., 0./ 

data dcp / 

4 . 208154e2, .580300e3, -.525250e3, -.642505e3, -.791376e3, 
4 -.182697e2, -.128857e2, -.422349e3, .176381e4, 

4 . 112521e2, . 160585e3, -.169299e2, -.532403e3, ~.238528e3, 
4 -.182071el, - . 850086el, -.295269e2, 0., 

4 . 365225e2, -.730849e3, -.502696e3, .370582e3, .517928el, 

4 .225676e2, -.216688e3, 0., 0., 

4 . 157475e2, -.627915e2, .757047e3, -.100614e3, .326629el, 

4 .595868el, .141182e3, -.682883e3, 0., 

4 .879873el, -.897230e2, .376565e3, -.403342e2, .314736el, 

4 -. 936668el, -.836769e2, 0., 0., 

4 . 591074el, .203259e2, .128924e3, -.102436e2, -.350311el, 

4 -.242638e2, .798459e3, 0., 0., 

4 . 519284el, .896935el, .381171e2, -.137585el, -.510951el, 

4 -.752038e2, 0., 0., 0./ 

data ecp / 

4 . 116592e2, .305728e3, -.241298e3, -.166628e3, -.102433e3, 
4 .138792el, .135712el, .482128e2, -.348874e3, 

4 . 570825el, .884544e2, .207647e2, -.119389e3, -.216169e2, 

4 .229104el, .253250el, -.792067el, 0., 



180195el, 


c 


c 


c 


c 


c 


c 


& . 203928e2, -.281133e3, 
& . 100733el, . 277132e2, 

& .820277el, -.154364e2, 
& . 821623, -.156018e2, . 
& . 412806el, -,28065162, 
& .987515, . 339060e2, 0. 

& .244269el, .100197e2, 

& .263659el, -.175990e3, 
& .202191el, . 384233el, 

& . 129598e2, 0., 0., 0./ 


- . 107793e3 
0 . , 0 . , 
.205793e3, 
143351e3, 0 
. 915792e2, 
, 0 ., 

.263694e2, 

0 . , 0 . , 

. 628850el, 


, . 457650e2, . 

168003e2, .838365, 
458728el, . 922570e-l 
-. 333185e-l, .184168, 
.743313, .877002, 


r 


data az / 

£ .710750, - . 614415el, -.632086e2, -.467833e2, .556705e2, 

& .824286, . 746758el, -.385889e2, -.455262e2, .809623e2, 

& .873086, - . 195828el, -.417508e2, -.431463e2, .208036e3, 

£ .904213, . 124751el, -.325326e2, -.428667e2, .217096e3, 

& . 102671el, . 387376e2, -.161621e2, -.255245e2, -.784807e2, 
& .970875, - . 10200el, -.993188el, .398457e-l, 0., 

& . 103304el, - . 555015el, .202955e2, 0., 0./ 


data bz / 

£ . 107229el, .861656el, .370722e2, .139011e2, -.135009e2, 

& .625098, - . 460729el, .209649e2, .121138e2, -.162146e2, 

& .434929, . 324383el, .199010e2, .101757e2, -.342626e2, 

£ .311295, .485004, .137742e2, .888031el, -.309522e2, 

& -.465274e-l, -.204439e2, .637080el, .419968el, .129796e2, 
& . 869030e-l, .186655el, .353080el, -.612253, 0., 

& -.585872e-l, .157079el, -.323532el, 0., 0./ 


data cz / 

& -.125673el, -.370256el, -.776456el, -.138693el, .118386el, 

& -.689867, . 109594el, -.398276el, -.108251el, .120105el, 

& -.454400, -. 123210el, -.334091el, -.804882, .210825el, 

& -.302086, -.321087, -.203163el, -.620696, .165245el, 

& . 972123e-2, .404607el, -.827695, -.208573, -.758996, 

& -.737745e-l, -.540958, -.389667, .997321e-l, 0., 

£ .237877e-l, -.115055, .203092, 0., 0./ 

data dz / 

£ .564944, .681208, .706484, .592861e-l, -.427210e-l, 

£ .286982, - . 898428e-l, .327436, .415356e-l, -.375039e-l, 

£ .176448, .198816, .243749, .274096e-l, -.563525e-l, 

£ .107468, . 632573e-l, .130377, .188157e-l, -.384201e-l, 

£ . 417402e-2, -.344141, .466769e-l, .395832e-2, .194343e-l, 

£ .218303e-l, .639254e-l, .187431e-l, -.411847e-2, 0., 

£ - .281715e-2, .324023e-2, -.525620e-2, 0., 0./ 

data ez / 

£ -. 822333e-l, -.443045e-l, -.233636e-l, -.903887e-3, .551468e-3, 
£ -.376727e-l, .162238e-2, -.970559e-2, -.569596e-3, .424122e-3, 

£ -.212727e-l, -.110471e-l, -.644569e-2, -.334336e-3, .555405e-3, 
£ - . 11692 4e-l, - . 364522e-2, -.302863e-2, -.206237e-3, .330019e-3, 
£ - . 536830e-3, .107287e-l, -.941988e-3, -.175392e-4, -.182292e-3, 
£ - . 17 97 62e-2, -.255478e-2, -.320898e-3, .542207e-4, 0., 

£ . 168221e-3, -.18832e-4, .489857e-4, 0., 0./ 

data amu / 

£ - . 1160076e-4, -.9105422, .1463029e-l, -.2140374e-2, 

£ .2397194e-4, -.5784272, .1658118e-l, .6903134e-2, 

£ . 5085043e-4, -.341487, .2450600e-l, - . 3561146e-l, 

£ . 6394112e-4, -.2376368, .6309492e-3, -.1622687el, 

£ . 5781887e-4, -.1844238, .2606784e-l, 0., 

£ .7256455e-4, - . 9524274e-l, .5037513e-l, 0., 

£ .7609039e-4, - . 7868582e-l, 0., 0./ 


c 


data bmu / 



fri fl-> fti (Mh(h(Mh(h(MHCi(h^fr'(h fii <h (H »Hh (h (MMH tfr #* th V> fr» IMMMU ^ (h (^ 


& . 6656010e-3, .4949794, - . 5019958e-2, .6529285e-3, 

& .5564725e-3, .2816531, - .5027652e-2, -. 1345295e-2, 

& .4774840e-3, .1473594, - . 6697224e-2, .7255623e-2, 

& . 4385020e-3, .9006170e-l, .6108099e-3, .3035173, 

& . 4438221e-3, .6040101e-l, - . 4562535e-2, 0., 

& . 4050530e-3, .2589951e-l, - . 8081647e-2, 0., 

& .3891948e-3, .1820922e-l, 0., 0./ 

data emu / 

-.2933969e-3, 

1970968e-3, 

-. 1322133e-3, 

1024141e-3, 

-. 1020840e-3, 

-.7626766e-4, 

-. 6458779e-4, 

data dmu / 

.7427050e-4, .1123425e-l, - . 4723839e-4, .3865996e-5, 

. 4272210e-4, .5058384e-2, - . 3715711e-4, -.4234384e-5, 

.2362256e-4, .2030404e-2, - . 4070960e-4, .2324839e-4, 

. 1654305e-4, .9370344e-3, .9381977e-5, ,8445985e-3, 

. 1688754e-4, .4609058e-3, - . 1018512e-4, 0., 

. 1114437e-4, .1227975e-3, - . 1682098e-4, 0., 

. 8791566e-5, .6257766e-4, 0., 0./ 

data emu / 

- . 6456605e-5, - . 5896774e-3, .1623374e-5, - . 9908122e-7, 
-.2690853e-5, - . 2352317e-3, .1135683e-5, .8514686e-7, 

-. 8014978e-6, - . 8216415e-4, .1134307e-5, - . 4590857e-6, 

- . 5014106e-6, - . 3279124e-4, -.2960969e-6, -.1570909e-4, 

- . 8622324e-6, - . 1377229e-4, .1576999e-6, 0., 

5020411e-6, - . 2772500e-5, .2731352e-6, 0., 

-.4216496e-6, - . 1234998e-5, 0., 0./ 

data fmu / 

. 8752161e-7, .1229026e-4, - .2239581e-7, .9916638e-9, 
-.3009241e-7, .4336410e-5, -. 1397984e-7, -.6893227e-9, 

-. 6458338e-7, .1314540e-5, -. 1276018e-7, .3600777e-8, 
-.3710875e-7, .4529650e-6, .3444222e-8, .1166667e-6, 
-.2239193e-9, .1623637e-6, - . 9011456e-9, 0., 

.7074486e-10, .2383398e-7, -.1794872e-8, 0., 

. 4509800e-8, .9579999e-8, 0., 0./ 

data ak / 

.395299el, .119879e2, -.832682e2, -.103603e4, .261125e2, 

.246095e2, -.571805e2, 

. 199665el, -.831120e2, -.110139e3, .299875e3, .434485e2, 

. 895136el, -.422029e2, 

. 198558el, .595832e2, -.442143e2, -.584437e3, .373716e2, 
-. 143675e2, .502985el, 

. 105928el, .101351e3, .830640el, -.318301e3, .469099e2, 

.154279e2, 0., 

.334316, . 109992e2, .124072e2, -.189644e3, .298795e2, 

. 844897el, 0., 

.413573, . 821184e2, .113875e2, -.723261e2, -.382696el, 

0., 0., 

.208749, . 378677e2, .223116e2, .792550el, 0., 0., 0./ 

data bk / 

. 386816e2, .412181e2, -.419438e3, -.242470e4, .411940el, 

-.507490e2, .168628e3, 

. 194822e2, -.560438e3, -.481050e3, .923042e3, .464790el, 
~.322183e2, .144838e3, 

& . 189164e2, .288748e3, -.206207e3, -.873106e3, -.115449e2, 
& . 801073el, - . 960227el, 


-.1060568, . 6886543e-3, -.7290226e-4, 

-. 5377449e-l, .6106363e-3, .1061916e-3, 

- . 2471167e-l, .7362709e-3, -.5837678e-3, 
-. 1315352e-l, - . 1286661e-3, -.2266401e-l, 
-.7566737e-2, .3111533e-3, 0., 

- . 2593217e-2, .5209350e-3, 0., 

-. 1543467e-2, 0., 0./ 



& . 10092 4e2, .490653e3, -.324274e2, -.306306e3, 

& -.541310e2, 0., 

& . 328202el, .387106e2, -.147438e2, -.828711e2, 

& ~.358117e2, 0., 

& .383393el, .308927e3, -.133907e2, .143656e2, . 
& 0., 0., 

& . 192122el f . 123284e3, .336369, -.216552e2, 0., 
data ck / 

& . 140687e3, .717156el, -.751764e3, -.206135e4, 

& . 369131e2, - . 181577e3, 

& . 706404e2, -.140314e4, -.757873e3, .992814e3, 

& .350726e2, -.182586e3, 

& . 668384e2, .500756e3, -.324643e3, -.445088e3, 

& . 146420e2, .196818el, 

& .356709e2, .868620e3, -.942568e2, -.782124e2, 

& . 693640e2, 0., 

& . 119939e2, .387282e2, -.530293e2, .998789el, . 
& .553921e2, 0., 

& . 131885e2, .423174e3, -.337860e2, .135247e2, - 
& 0., 0., 

& . 658813el, .144224e3, -.142705e2, .204578e2, 0 
data dk / 

& . 226110e3, -.911924e2, -.566912e3, -.757541e3, 
& 897288el, .859388e2, 

& . 113538e3, -.154128e4, -.505860e3, .442621e3, 

& -.129798e2, .101698e3, 

& . 104546e3, .363789e3, -.204871e3, -.927269e2, 

& ~.117248e2, .643952el, 

& .561818e2, .666792e3, -.647282e2, -.466313el, 

& ~.366810e2, 0., 

& .200944e2, .548304el, -.299886e2, .227739el, . 
& -.353787e2, 0., 

& .207305e2, .250668e3, -.122339e2, -.233991el, 

& 0 . , 0 * , 

& . 107630e2, .728083e2, -.134534el, -.597164el, 
data ek / 

& . 127138e3, -.810415e2, -.157470e3, -.108281e3, 
& 623025el, -.213335e2, 

& .599079e2, -.632398e3, -.125800e3, .634709e2, 

& -.498154el, -.270417e2, 

& . 527822e2, .844428e2, -.494556e2, -.128529e2, 

& -.389761el, -.896353el, 

& .249670e2, .180596e3, -.180857e2, -.585083el, 

& . 115271el, 0., 

& . 462882el, -.120106e2, -.961485el, -.581069el, 
& .274595el, 0., 

& . 427728el, .475889e2, -.610064el, -.556444el, 

& 0., 0., 

& -.127699el, .684807el, -.498832el, -.485454el, 
data apr / 

& .7695318, . 4081926e2, -.5060907e3, .3172744e2, 
& . 1034825e3, -.3414530e3, .1398841e2, 

& .7483071, . 4355608e2, .6877578e3, -.7024989e3, 
& . 9153227e2, .2304245e4, -.1424413e2, 

& .7406012, . 3206825e2, .4409105e3, .2474049e2, 

& . 9436717e2, -.7123564e2, -.7846718e2, 

& .7548377, .6190687e2, -.2641868e3, .7451996e2, 
& -.1252168e3, -.5022975e4, 0., 

& .8037721, . 9173178e2, .1874806e2, -.1974460e3, 
& . 1619271e3, -.1396134e3, 0., 

& .8167207, - . 1857130e2, -.4287169e2, .6447596e2 
& 0., 0., 0., 


-.330961e2, 

-.381078e2, 

146502e2, 

0 ., 0 ./ 

-. 186054e2, 
-. 155778e2, 
113653e2, 
-. 146607el, 
117041e2, 

. 187337e2, 

- , 0 ., 0 ./ 

-. 645054el, 
-.220224el, 

. 135799el, 

. 306898el, 
122011el, 
.107119e2, 
0 ., 0 ., 0 ./ 

-. 621476el, 
-. 558790el, 
-.542822el, 
-.562490el, 
-.578171el, 
-.717162el, 
0 ., 0 ., 0 ./ 

-.7624622e4 
-.3107328e3 
-. 9639789e3, 
. 1034771e3, 
.5713095e2, 
, .7384318e2 



c 


& .814777, - . 6650el, -.2736096e2, .4889586e2, 

£ 0 ., 0 ., 0 ./ 


5525036e2, 


c 


c 


c 


data bpr / 

£ -.2487156, -. 6945482e2, .5765441e3, -.1923593e2 
£ -.3861297e2, .1193445e3, -.2209182el, 

& -.1659050, -. 6751147e2, -.5984367e3, .4515788e3 
& -.3041980e2, -.6078459e3, .3715427el, 

£ -.1629686, - . 4531855e2, -.3447078e3, .2323747el 
£ -.2761116e2, .2178196e2, .1520447e2, 

£ -.2611566, - . 8052889e2, .2214640e3, -.2929220e2 
£ .2959573e2, .9203587e3, 0., 

£ -.5531739, - . 1294114e3, -.1571672e2, .9815976e2 
£ -.3657252e2, .2786187e2, 0., 

£ -.618678, .2610712e2, .3082544e2, -.2113068e2, 

£ 0 . , 0 . , 0 . , 

£ -.6129922, . 9087650el, .1729215e2, -.1351549e2, 

£ 0 ., 0 ., 0 ./ 


, .3908646e4, 
, . 1863582e3, 
, .3960004e3, 
, -.2138868e2 
, -,2693925e2 
-. 1345184e2, 

. 1123889e2, 


t 


t 


data cpr / 

£ -.3371200e-2, .4642510e2, -.2594581e3, .6031312el, -.7963915e3, 

£ . 5742599el, -.1639670e2, .1254939, 

£ -.3548542e-l, .4101816e2, .2059622e3, -.1153420e3, -.3984823e2, 
£ . 4035104el, .6397478e2, -.3586053, 

£ . 8014385e-l, .2519178e2, .1064873e3, -.4089151el, -.6388451e2, 

£ .3229296el, -.2554952el, -.1165686el, 

£ .3168480, . 4117452e2, -.7330522e2, .4227110el, .1519602el, 

£ -.2775372el, -.6737887e2, 0., 

£ .9090234, . 7264736e2, .5656501el, -.1915726e2, .4625149el, 

£ .3311637el, -.2213081el, 0., 

£ . 1015263el, -.1366805e2, -.8483934el, .2725852el, .9862360, 

£ 0 ., 0 ., 0 ., 

£ . 1017670el, - . 4283159el, -.4132288el, .1493734el, -.8807574, 

£ 0 ., 0 ., 0 ./ 


data dpr / 

£ .5237761, - . 1500926e2, .5768434e2, -.1191714el, .8069517e2, 

£ -.4251452, . 1109347el, -.2782377e-2, 

£ .3843113, -. 1204527e2, -.3502408e2, .1460894e2, .4000552el, 

£ -.2667702, - . 3358410el, .1651461e-l, 

£ .1480977, - . 6755087el, -.1623127e2, .8375614, .5082427el, 

£ -.1885223, .1451363, .4430242e-l, 

£ -.1019887, - . 1025306e2, .1198777e2, -.2590555, -.3418909e-l, 

£ .1290359, .2463703el, 0., 

£ -.6243393, -.2015462e2, -.1040290el, .1841274el, -.3686820, 

£ -.1501584, .8742029e-l, 0., 

£ -.6989158, . 3448930el, .1132069el, -.1705427, - . 3633532e-l, 

£ 0 . , 0 . , 0 . , 

£ -.7111724, .9588901, .4767165, - . 8122753e-l, .3371865e-l, 

£ 0., 0., 0./ 

data epr / 

£ -.4272376, .2350767el, -.6333573el, .1301576, -.4069244el, 

£ . 1565525e-l, -. 3701973e-l, .6736752e-5, 

£ -.2847560, . 1712608el, .2943975el, -.9148849, -.1930037, 

£ . 8781471e-2, .8795057e-l, - . 3685522e-3, 

£ -.1294933, .8752933, .1221500el, - . 6543669e-l, -.2000122, 

£ .5488061e-2, - . 4016437e-2, - . 8360113e-3, 

£ -.1326718e-l, .1246034el, -.9677752, .5423023e-2, -.5284942e-3, 

£ -.2973274e-2, - . 4499464e-l, 0., 

£ .1903222, .2763213el, .9598832e-l, - . 8709938e-l, .1391085e-l, 

£ .3407133e-2, -. 1716422e-2, 0., 

£ .2151091, -.4205156, - . 7326105e-l, .5177594e-2, .6720911e-3, 

£ 0 . , 0., 0., 

£ .2215290, -.1026406, -.2655841e-l, .2173037e-2, - . 6355462e-3, 

£ 0., 0., 0./ 



data fpr / 

& .9284843e-l, -.1430131, .2747920, - . 5625333e-2, .8174707e-l, 

& -.2291528e-3, .4881128e-3, .3740205e-6, 

& . 5734562e-l, -.94534676-1, - . 9787067e-l, .2261774e-l, . 3624133e-2, 
& - . 1150000e-3, - . 9193691e-3, .3216955e-5, 

& .2544039e-l, - . 4395200e-l, - . 3631877e-l, .1793128e-2, .3121600e-2, 

& - . 6365539e-4, .4349632e-4, .6273942e-5, 

& . 6541023e-2, - . 592 6933e-l, .3085867e-l, .1836850e-4, .2331795e-4, 

& .2716070e-4, .3283492e-3, 0., 

& -.2182155e-l,-. 1497413,-. 3499282e-2, . 1622326e-2, -.201543 6e-3, 

& -.3092820e-4, .1339615e-4, 0., 

& - .2479595e-l, .1987815e-l, . 1844057e-2, - . 6127611e-4, - . 4987500e-5, 

& 0., 0. , 0., 

& - .2557178e-l, .4226421e-2, . 5732041e-3, - . 2296542e-4, .4740705e-5, 

& 0., 0., 0./ 
c 

c if temperature is less than 500K 
c 

if (tinp .It. 500. .or. iflag .eq. 0) then 
cp = .24 
z = 1 . 

xmu * 1 . 4584e-5*tinp**l . 5/ (tinp+110 . 33) 
xk = 5 . 9776e-6*tinp**l . 5/ (tinp+194 . 4) 
pr = .24* xmu / xk 

rho * pinp*airmO/ (runiv*tinp) *2 . 4218e-2 
airm = airmO 

gamma = 1 . / (1 .-runiv/airm/cp) 
return 
endif 
c 

c Find pressure range for interpolation 
c 

if (pinp .It. l.e-4 .and. pinp .gt. 100.) then 

write < 6,*) ' Input pressure of ',pinp,' is outside range of', 

& 9 available curve fit' 

ifatal * 1 
return 
endif 
c 

c Find pressures which bound the input pressure 
c 

do 10 i = 1,7 
if (pinp .le. ptb(i)) then 
if (pinp .eq. ptb(i)) then 
intp = 0 
else 

intp * 1 
endif 
ip * i 

ipml ~ ip -1 
go to 20 
endif 

10 continue 
c 

write (6,*) ' Pressure is outside the range of available data' 

ifatal * 1 

return 

20 continue 
c 

xi = alog (tinp* . 0001) 
c 

c Specific heat 
c 

do 30 i = l,ncpmax(ip) 

if (tinp .gt. tcpmax (i, ip) ) go to 30 



it * i 
go to 40 
30 continue 
c 

write (6,*) ' Temperature ',tinp, ' is outside the range of ' , 
available data - cp high' 
ifatal = 1 
return 

40 continue 

cpm2 = acp (it, ip) *xi**4 + bcp (it , ip) *xi**3 + ccp(it,ip)* 

& xi**2 + dcp (it, ip) *xi + ecp(it,ip) 

c 

if(intp .eq. 0 .or. ipml . le. 0) then 
cp = exp ( cpm2 ) 
go to 65 
endif 

do 50 i * 1 , ncpmax ( ipml ) 
if (tinp .gt. tcpmax (i, ipml ) ) go to 50 
it ® i 
go to 60 
50 continue 
c 

write (6,*) 9 Temperature ',tinp, 'is outside the range of', 

&' available data - cp low' 
ifatal = 1 
return 

60 continue 

cpml = acp (it, ipml) *xi**4+bcp (it, ipml) *xi**3+ccp (it , ipml) *xi**2 
& + dcp (it, ipml) *xi + ecp(it,ipml) 

c 

c interpolation for specific heat 
c 

alogcp = (cpm2 - cpml) / (alog (ptb (ip) ) -alog (ptb (ipml) )) * 

& (alog(pinp) - alog (ptb (ipml) ) ) + cpml 

cp * exp (alogcp) 
c 

c Compressibility factor 
c 

65 continue 

xi = tinp*. 001 
do 70 i = l,nzmax(ip) 
if (tinp .gt. tzmax(i,ip)) go to 70 
it * i 
go to 80 
70 continue 
c 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - z high' 
ifatal = 1 
return 

80 continue 

z2 = az(it,ip) + bz(it,ip)*xi + cz (it, ip) *xi**2 
& + dz (it, ip) *xi**3 + ez (it , ip) *xi**4 

if(intp .eq. 0 .or. ipml .le. 0) then 
z = z2 
go to 105 
endif 

c 

do 90 i = 1, nzmax (ipml) 
if (tinp .gt. tzmax (i, ipml) ) go to 90 
it = i 
go to 100 
continue 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - z low' 


90 



ifatal « 1 
return 

100 continue 

zl « az (it, ipml) + bz (it, ipml) *xi + cz (it , ipml) *xi**2 
& + dz (it , ipml) *xi**3 + ez (it , ipml) *xi**4 

c 

c interpolation of compressibility factor 
c 

z * (z2 - zl)/(ptb(ip) - ptb(ipml)) * (pinp - ptb(ipml)) + zl 
c 

c Viscosity 
c 

105 continue 

do 110 i * l,nmumax(ip) 
if (tinp .gt. tmumax (i, ip) ) go to 110 
it * i 
go to 120 
110 continue 
c 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - mu high' 
ifatal = 1 
return 

120 continue 

xmu2 * amu(it,ip) + bmu (it, ip) *xi + emu (it , ip) *xi**2 
& + dmu (it, ip) *xi**3 + emu (it, ip) *xi** 4 + fmu (it, ip) *xi**5 

if(intp .eq. 0 .or. ipml .le. 0) then 
xmu = xmu2 
go to 145 
endif 
c 

do 130 i = 1 , nmumax ( ipml ) 
if (tinp .gt. tmumax ( i , ipml ) ) go to 130 
it * i 
go to 140 
130 continue 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - mu low' 
ifatal * 1 
return 

140 continue 

xmul = amu(it,ipml) + bmu (it, ipml) *xi + emu (it, ipml) *xi**2 
& + dmu (it, ipml) *xi**3 + emu (it, ipml) *xi**4 + fmu (it, ipml) *xi**5 

c 

c interpolation of viscosity 
c 

xmu = (xmu2 - xmul) / (ptb (ip) - ptb(ipml)) * (pinp - ptb(ipml)) 

& + xmul 

c 

c Thermal conductivity 
c 

145 continue 

xi * alog (tinp* . 0001) 

do 150 i = l,nkmax(ip) 
if (tinp .gt. t3onax (i, ip) ) go to 150 
it * i 
go to 160 
150 continue 

write (6,*) ' Temperature ',tinp,' is outside the range of', 

&' available data - K high' 
ifatal = 1 
return 

160 continue 

xk2 = ak (it , ip) *xi**4 + bk (it, ip) *xi**3 + ck(it,ip)* 


c 



& xi**2 + dk (it, ip) *xi + ek(it,ip) 

if(intp .eq. 0 .or. ipml . le. 0) then 
xk = exp <xk2 ) 
go to 185 
endif 
c 

do 170 i = 1 , nkmax ( ipml ) 
if (tinp .gt. tkmax (i, ipml) ) go to 170 
it - i 
go to 180 
170 continue 
c 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - K low' 
ifatal = 1 
return 

180 continue 

xkl = ak (it, ipml) *xi**4 + bk (it , ipml) *xi**3 + ck (it, ipml) *xi**2 
& + dk (it, ipml) *xi + ek (it, ipml) 

c 

c interpolation for thermal conductivity 
c 

alogk = (xk2 - xkl) / (alog (ptb (ip) ) -alog (ptb (ipml) )) * 

& (alog(pinp) - alog (ptb (ipml) ) ) + xkl 

xk = exp (alogk) 
c 

c Prandtl number 
c 

185 continue 

xi = tinp *.001 
do 190 i = l,nprmax(ip) 
if (tinp .gt. tprmax (i, ip) ) go to 190 
it = i 
go to 200 
190 continue 
c 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - Pr high' 
ifatal = 1 
return 

200 continue 

pr2 = apr ( it , ip) + bpr (it , ip) *xi + cpr (it , ip) *xi**2 
& + dpr (it , ip) *xi**3 + epr (it , ip) *xi**4 + fpr (it , ip) *xi**5 

if (intp .eq. 0 .or. ipml . le . 0) then 
pr = pr2 
go to 300 
endif 
c 

do 210 i = l,nprmax(ipml) 
if (tinp .gt. tprmax (i, ipml) ) go to 210 
it = i 
go to 220 
210 continue 

write (6,*) ' Temperature ',tinp, 'is outside the range of', 

&' available data - Pr low' 
ifatal = 1 
return 

220 continue 

prl = apr (it, ipml) + bpr (it , ipml ) *xi + cpr (it , ipml) *xi**2 
& + dpr (it, ipml) *xi** 3 + epr ( it , ipml ) *xi**4 + fpr (it , ipml) *xi**5 

c 

c interpolation of Prandtl number 
c 

pr = (pr2 - prl) / (ptb (ip) - ptb(ipml)) * (pinp - ptb(ipml)) 

& + prl 

c 



300 continue 
c 

c compute density 
c 

rhoO = pinp*airm0/ (runiv*tinp) *2 . 4218e-2 
if (z . It . 1 . ) z = 1 . 
rho « rhoO/z 
c 

c compute molecular weight and effective gamma 
c 

airm - airmO/z 

gamma - 1 . / (1 . -runiv/airm/cp) 
c 

write (6,*) 9 cp = ',cp 
write (6, *) 9 z - 9 , z 
write (6,*) 9 rho = ',rho 
write (6,*) ' mu = ',xmu 
write (6, *) ' k = ' ,xk 
write (6,*) 9 pr = ',pr 

c 

return 

end 

c 

subroutine find(x0, yO, zO) 
c 

c Given a location (x0,y0,zQ), and the grid cell containing this point, 
c this routine uses Newton's method to find corresponding location 
c (xi,eta, zeta) in the master element, 
c 

common /pos/ xi, eta, zeta 

common /psidps/ psi(8), dpdxi(8), dpdeta(8), dpdzta(8) 
common /connty/ il,jl,kl 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151,40,110), blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

common /save/ igsav, issav,isav, jsav,ksav,asav,bsav, gsav, 

& first, back, second, ifatal 

logical blank, first, second, back 
dimension x(8), y(8), z(8), r(3), fj(3,3), dr (3) 
c 

c set initial guess at xi=0, eta^O and zeta=0 
c 

write(6,*) 'x0,y0,z0 = ',x0,y0,z0 

xi = 0. 

eta = 0 . 

zeta * 0 . 

ipl = il + 1 

jpi - ji + i 

kpl - kl + 1 
c 

x(l) = xyz (il, jl,kl, 1) 
y (1) = xyz (il, jl,kl,2) 
z(l) = xyz (il, jl,kl, 3) 
c 

x (2) = xyz (ipl, jl,kl, 1) 
y (2) = xyz (ipl, jl,kl,2) 
z (2) = xyz (ipl, jl,kl,3) 
c 

x(3) = xyz (ipl, jl, kpl, 1) 
y (3) = xyz (ipl, jl, kpl, 2) 
z(3) = xyz (ipl, jl, kpl, 3) 
c 

x (4) = xyz (il, jl/kpl, 1) 
y (4) = xyz (il, jl,kpl,2) 
z (4) = xyz (il, jl,kpl, 3) 
c 



x<5) - xyz (il, jpl,kl, 1) 
y (5) - xyz (il, jpl,kl,2) 
z (5) = xyz (il, jpl,kl,3) 

x(6) = xyz (ipl, jpl,kl, 1) 
y(6) = xyz (ipl, jpl, kl, 2) 
z (6) = xyz (ipl, jpl, kl, 3) 

x (7 ) = xyz (ipl, jpl,kpl, 1) 
y (7) = xyz (ipl, jpl, kpl, 2) 
z (7) * xyz (ipl, jpl, kpl, 3) 

x (8) = xyz (il, jpl, kpl, 1) 
y (8) = xyz (il, jpl, kpl, 2) 
z (8) = xyz (il, jpl, kpl, 3) 

iter - 0 
10 continue 


c 

c compute shape function and it's derivatives 


c 

call shape 


c 

c initialize jacobian matrix 
c 

do 20 j = 1,3 
do 20 i=l,3 
f j (if j) = 0 . 

20 continue 


r (1) = xO 
r (2) = yO 
r(3) = zO 


c 

c get jacobian matrix 


c 


c 


do 30 i 
f j(l,l) 
f j (2, 1) 
f j (3, 1) 
f j (1,2) 
f j(2,2) 
f j (3,2) 
f j (1,3) 
f j (2,3) 
f j(3,3) 
r (1) = 
r (2) = 
r (3) = 


- 1,8 
= fj(l,D 
= f j (2,1) 

= f j (3, 1) 

= f j (1,2) 

= f j (2,2) 

= f j (3,2) 

= f j (1, 3) 

= f j (2, 3) 

= f j (3, 3) 

(1) - x (i) 

(2) - y(i) 


x (i 

y(i 

z (i 
x (i 

y(i 

z (i 
x(i 

y(i 

z (i 
*psi ( 
*psi ( 
*psi ( 


) *dpdxi (i) 

) *dpdxi (i) 

) *dpdxi (i) 

) *dpdeta (i) 
) *dpdeta (i) 
) *dpdeta (i) 
) *dpdzta (i) 
) *dpdzta (i) 
) *dpdzta (i) 

i) 

i) 

i) 


30 continue 

det = f j(l,l)*fj(2,2)*f j (3, 3) + f j (1,2) *f j (2,3) *f j (3, 1) 
& + f j (1,3) *f j (2,1) *f j (3,2) - f j(3,l)*f j(2,2)*f j(l,3) 

& - f j (1, 2) *f j(2,l)*f j (3, 3) - f j(l,l)*f j(3,2)*f j(2,3) 

write (6,*) ' determinant of Jacobian = ',det 
res = sqrt(r(l) *r(l) + r(2)*r(2) + r(3)*r(3)) 
if (res .It. 1.0e-4) go to 9999 
idm = 3 


call gauss (fj, r, dr, idm) 
if (ifatal .ne. 0) return 
c write (6,*) dr 

sdr = abs(dr(l)) + abs(dr(2)) + abs(dr(3)) 
xi = xi + dr(l) 
eta = eta + dr (2) 
zeta = zeta + dr (3) 
if (sdr ,ge. l.e-3) then 
iter = iter + 1 



if (iter .gt. 40) then 

write (6,*) 9 !!! iteration limit exceeded !!!' 

write (6,*) 9 sdr = f ,sdr 
ifatal - 1 
return 
else 

go to 10 
endif 
endif 

9999 continue 

if(abs(xi) .gt. 1.) xi - sign(l.,xi) 
if(abs(eta) .gt. 1.) eta = sign {!., eta) 
if(abs(zeta) .gt. 1.) zeta = sign (1 . , zeta) 
call shape 

write (6,*) ' xi, eta, zeta = ' , xi, eta, zeta 
return 
end 
c 

subroutine shape 
c 

c This routine evaluate 3-D linear shape functions and their 
c derivatives, 
c 

common /pos/ xi, eta, zeta 

common /psidps/ psi(8), dpdxi(8), dpdeta(8), dpdzta(8) 
c 

omxi * 1 . - xi 
opxi = 1 . + xi 
ometa = 1. - eta 
opeta = 1. + eta 
omzeta * 1 . - zeta 
opzeta * 1 . + zeta 
c 

omxiet = omxi * ometa 

opxime * opxi * ometa 

opxiet * opxi * opeta 

omxipe * omxi * opeta 

c 

psi(l) = .125 * omxiet * omzeta 

psi(2) = .125 * opxime * omzeta 

psi(3) * .125 * opxiet * omzeta 

psi(4) = .125 * omxipe * omzeta 

psi(5) * .125 * omxiet * opzeta 

psi<6) = .125 * opxime * opzeta 

psi(7) - .125 * opxiet * opzeta 

psi(8) * .125 * omxipe * opzeta 

c 

ometze = ometa * omzeta 
ometpz * ometa * opzeta 
opetze = opeta * opzeta 
opetmz = opeta * omzeta 
omxize = omxi * omzeta 

omxipz * omxi * opzeta 

opxize « opxi * opzeta 

opximz = opxi * omzeta 

c 

dpdxi(l) = -.125 * ometze 
dpdxi (2) = -dpdxi(l) 
dpdxi(3) = .125 * opetmz 
dpdxi { 4 ) = -dpdxi ( 3 ) 
dpdxi (5) = -.125 * ometpz 
dpdxi (6) * -dpdxi (5) 
dpdxi (7) = .125 * opetze 
dpdxi (8) = -dpdxi (7) 
c 

dpdeta(l) = -.125 * omxize 



dpdeta(2) = -.125 * opximz 
dpdeta(3) * -dpdeta(2) 
dpdeta(4) - -dpdeta(l) 
dpdeta(5) * -.125 * omxipz 
dpdeta(6) * -.125 * opxize 
dpdeta(7) * -dpdeta(6) 
dpdeta(8) * -dpdeta(5) 
c 

dpdzta(l) = -.125 * omxiet 
dpdzta(2) = -.125 * opxime 
dpdzta(3) * -.125 * opxiet 
dpdzta(4) = -.125 * omxipe 
dpdzta(5) = -dpdzta(l) 
dpdzta(6) = -dpdzta(2) 
dpdzta(7) = -dpdzta(3) 
dpdzta(8) - -dpdzta(4) 
c 

return 

end 

c 

SUBROUTINE GAUSS (F JAC, C, X, IDIM) 

common /save/ igsav, issav, isav, jsav, ksav, asav,bsav, gsav, 

& first, back, second, ifatal 

C 

C GAUSS ELIMNATION SOLVER FOR SYSTEM AX=C 
C 

c IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION F JAC (3,3) ,c(3) ,X(3) 

C 

C FORWARD ELIMINATION 
C 

DO 30 I = 1, IDIM-1 
IP1 =1+1 

IF (FJAC (1,1) .EQ. 0.) THEN 
WRITE(6, 9006) I 
ifatal = 1 
return 
ENDIF 

DO 20 J = IPl, IDIM 

IF (FJAC (J, I) .EQ. 0.) GO TO 20 

DO 10 K = IPl, IDIM 

FJAC ( J, K) = F JAC ( J, K) - FJAC (I, K) *FJAC ( J, I) /FJAC (I, I) 

10 CONTINUE 

C(J) = C(J) -C(I) *FJAC(J,I) /FJAC (I, I) 

20 CONTINUE 
30 CONTINUE 
C 

C BACK SUBSTITUTION 
C 

DO 70 I = IDIM, 1,-1 
X(I) = 0. 

DO 60 J = I, IDIM 

C(I) = C(I) - FJAC ( I, J) *X ( J) 

60 CONTINUE 

X(I) * C(I) / FJAC (1,1) 

70 CONTINUE 
C 

9006 FORMAT ( 9 *** ZERO PIVOT, I « ',13) 

RETURN 

END 

SUBROUTINE CELL3 ( IS, IE, JS, JE, KS, KE, 

C I, J, K, A, B, G, X, Y, Z , AMAT, BMAT, STATUS) 

************************************************************************ 

C 

C Find the point (x,y,z) in the grid XYZ and return its (i, j,k) cell 




non nnnnnn no o o ooooooo 


number. We ASSUME that the given (i,j,k), (a,b,g), and subset are 

valid. These can be checked with STRTxx. 

STATUS-1 - Unable to find the point without going out of the 

computational domain or active subet. The computat- 
ional point returned indicates the direction to look.. 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151, 40,110), blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

DIMENSION AMAT (3,3) , BMAT (3,3) 

LOGICAL BLANK 
INTEGER STATUS 

PARAMETER (MITER=5) 

LOGICAL DIALOG 
INTEGER DI,DJ,DK 

DATA I S AVI , JS AVI , KS AVI / 3 * 0 / , ISAV2 , JSAV2 , KS AV2 / 3 * 0 / 

DATA DIALOG/ .FALSE./ 

DATA DIALOG/. true./ 

STATUS- 0 

Reset saved points used for checking if we're stuck on the same 
point (or ping-ponging back and forth) when searching for the right 
cell. This would happen while following a point outside the 
computational domain. 

ISAV1- 0 
ISAV2- 0 
DI= 1 
DJ= 1 
DK= 1 

Maximum number of steps before we'll give up is max (IDIM, JDIM, KDIM) . 

MSTEP- MAX (IDIM, JDIM, KDIM) 

NSTEP— 0 

20 CONTINUE 

NSTEP— NSTEP+1 
C 

if(i .ge. idim) i = idim - 1 

if ( j . ge . jdim) j — jdim - 1 

if (k .ge. kdim) k = kdim - 1 

II = 1+1 
J1 = J+l 
K1 = K+l 
C 

XI = XYZ (I, J, K, 1) 

X2 = XYZ (II, J,K, 1) 

X3 = XYZ (I, J1,K, 1) 

X4 = XYZ (II, J1,K,1) 

X5 = XYZ (I, J, Kl, 1) 

X6 = XYZ (II, J,K1,1) 

X7 = XYZ (I, Jl, Kl, 1) 

X8 = XYZ (II, J1,K1, 1) 

C 

Y1 = XYZ (I, J, K, 2 ) 

Y2 = XYZ (II, J,K,2) 

Y3 = XYZ ( I , Jl , K, 2) 

Y4 = XYZ (I1,J1,K, 2) 

Y5 = XYZ (I, J,K1,2) 

Y6 = XYZ (II, J, Kl, 2) 

Y7 = XYZ (I,J1,K1,2) 



c 


c 


c 


c 


c 


c 


c 


c 


Y8 = XYZ(I1, J1,K1,2) 

Z1 = XYZ (I, J, K, 3) 

Z2 = XYZ(I1, J,K,3) 

Z3 = XYZ (I, Jl, K r 3) 

Z4 - XYZ (Hr Jl, K, 3) 

Z5 - XYZ(I, J,K1,3) 

Z6 - XYZ (II, J,K1,3) 

Z7 - XYZ (I, Jl, Kl, 3) 

Z8 - XYZ (II, Jl, Kl, 3) 

XO - XI 

XA - X2-X1 

XB * X3-X1 

XG - X5-X1 

XAB •= X4-X3-X2+X1 

XAG = X6-X5-X2+X1 

XBG - X7-X5-X3+X1 

XABG= X8-X7-X6+X5-X4+X3+X2-X1 

YO = Y1 

YA = Y2-Y1 

YB = Y3-Y1 

YG = Y5-Y1 

YAB - Y4-Y3— Y2+Y1 

YAG = Y6-Y5-Y2+Y1 

YBG = Y7-Y5— Y3+Y1 

YABG- Y8-Y7— Y6+Y5-Y4+Y3+Y2-Y1 

ZO = Z1 

ZA - Z2-Z1 

ZB - Z3-Z1 

ZG = Z5-Z1 

ZAB - Z4-Z3-Z2+Z1 

ZAG - Z6-Z5-Z2+Z1 

ZBG = Z7-Z5-Z3+Z1 

ZABG= Z8-Z7-Z6+Z5-Z4+Z3+Z2-Z1 


A - .5 

B = .5 

G = .5 

ITER- 0 

30 CONTINUE 

ITER- ITER+1 

XH - XO+XA*A+XB*B+XG*G+XAB*(A*B)+XAG*(A*G)+XBG*(B*G) 

C +XABG* (A* (B*G) ) 

YH - Y 0 +YA* A+YB *B+YG *G+YAB* (A*B) +YAG* (A*G) +YBG* (B*G) 

C +YABG* (A* (B*G) ) 

ZH = ZO+ZA*A+ZB*B+ZG*G+ZAB* (A*B) +ZAG* (A*G) +ZBG* (B*G) 

C +ZABG* (A* (B*G) ) 


AMAT (1,1)= XA + XAB*B + XAG*G + XABG* (B*G) 
AMAT ( 2 , 1 ) = YA + YAB*B + YAG*G + YABG* (B*G) 
AMAT (3,1)= ZA + ZAB*B + ZAG*G + ZABG* (B*G) 


AMAT (1,2)= XB + XAB*A + XBG*G + XABG* (A*G) 
AMAT (2,2)= YB + YAB*A + YBG*G + YABG* (A*G) 
AMAT (3,2)= ZB + ZAB* A + ZBG*G + ZABG* (A*G) 


AMAT (1,3)- XG + XAG* A + XBG*B + XABG* (A*B) 
AMAT (2,3)= YG + YAG*A + YBG*B + YABG* (A*B) 
AMAT (3,3)= ZG + ZAG* A + ZBG*B + ZABG* (A*B) 


CALL INV3X3 (AMAT, BMAT, ISTAT) 

IF ( ISTAT. NE.O) THEN 

IF (DIALOG) WRITE (6,*) 'DEGENERATE VOLUME AT INDEX ',I,J,K 
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See if we're at the edge of the cell. If so, move away from the 
(possibly degenerate) edge and recompute the matrix. 


AS = A 
BS = B 
GS - G 

IF (A.EQ.O.) A = .01 

IF (A.EQ.l.) A - .99 

IF (B.EQ.O.) B « .01 

IF (B.EQ.l.) B - .99 

IF (G.EQ.O.) G - .01 

IF (G.EQ.l.) G - .99 

IF (A.NE.AS .OR. B.NE.BS .OR. G.NE.GS) THEN 
GOTO 30 


We're inside a cell and the transformation matrix is singular. Move 
to the next cell and try again. 


ELSE 

A = DI+.5 
B = DJ+.5 
G = DK+ . 5 
GOTO 40 
ENDIF 
ENDIF 


DX 

DY 

DZ 

DA 

DB 

DG 

A 

B 

G 


X-XH 

Y-YH 

Z-ZH 

DX*BMAT (1, 1) 

DX*BMAT (2,1) 

DX*BMAT (3, 1) 

A+DA 

B+DB 

G+DG 


+ DY*BMAT (1,2) 
+ DY*BMAT (2,2) 
+ DY*BMAT (3,2) 


+ DZ*BMAT (1, 3) 
+ DZ*BMAT (2, 3) 
+ DZ*BMAT (3, 3) 


If we're WAY off, don't bother with the error test. In fact, go 
ahead and try another cell. 


IF (ABS(A— .5) .GT.3. .OR. ABS (B- . 5) .GT . 3 . 

C .OR. ABS (G-. 5) .GT.3.) THEN 

GOTO 40 


Check iteration error and branch out if it's small enough. 


ELSE 

ERR2= DA*DA + DB*DB + DG*DG 
IF (ERR2/3.LE.1.E-4) GOTO 40 
ENDIF 

IF ( ITER. LT. MITER) GOTO 30 
40 CONTINUE 


The point is in this cell . 

IF (ABS (A-.5) .LE. .50005 .AND. ABS (B- . 5) . LE . . 50005 
C .AND. ABS(G-.5) .LE. .50005) THEN 

IF (DIALOG) WRITE (6,*) 'MATCH' 

We've taken more steps then we're willing to wait... 

ELSE IF (NSTEP .GT.MSTEP) THEN 
STATUS= 1 

IF (DIALOG) WRITE (6,*) 'MORE THAN ',MSTEP,' STEPS' 



C Update our (i,j,k) guess, keeping it inbounds. 

C 

ELSE 

IN = I 
JN = J 
KN = K 

IF (A.LT.O.) IN - MAX < IN-1, IS ) 

IF (A.GT.l.) IN - MIN(IN+1, IE-1) 

IF (B.LT.O.) JN - MAX (JN-1, JS ) 

IF (B.GT.l.) JN = MIN ( JN+1, JE-1) 

IF (G.LT.O.) KN - MAX (KN-1, KS ) 

IF (G.GT.l.) KN - MIN (KN+1,KE-1) 
c IF (g.LT.O.) JN - MAX (JN-1, JS ) 

C IF (g.GT.l.) JN - MIN (JN+1, JE-1) 

c IF (b.LT.O.) KN - min (KN+1,KS ) 

c IF (b.GT.l.) KN - max (KN-1, KE-1) 

IF (DIALOG) WRITE (*,*) 'TRY CELL INDEX ',IN,JN,KN 
C 

C Check IBLANK for this cell. 

C 

IF (BLANK) THEN 
INI = IN+1 

JN1 - JN+1 

KN1 = KN+1 

IF ( IBLANK ( IN, JN,KN) .EQ.O .OR. IBLANK ( INI, JN, KN) .EQ.O 

C .OR. IBLANK (IN, JN1,KN) .EQ.O .OR. IBLANK (INI, JN1,KN) .EQ.O 

C .OR. IBLANK ( IN, JN,KN1) .EQ.O .OR. IBLANK (INI, JN, KN1) .EQ.O 

C .OR. IBLANK (IN, JN1 , KN1 ) . EQ . 0 .OR. IBLANK (INI, JN1, KN1) . EQ . 0) 

C THEN 

STATUS- 1 

IF (DIALOG) WRITE (6,*) 'EXTRAPOLATE' 

GOTO 50 
END IF 
END IF 
C 

C Not repeating a previous point. Use the new (i,j,k) and try again. 

C 

IF ( (IN.NE. ISAV1 .OR. JN.NE.JSAV1 .OR. KN.NE.KSAV1) 

C .AND. (IN.NE. ISAV2 .OR. JN.NE.JSAV2 .OR. KN . NE . KS AV2 ) ) THEN 

ISAV2- IS AVI 
JSAV2= JSAV1 
KSAV2- KSAV1 
ISAV1= IN 
JSAV1= JN 
KSAV1= KN 

DI = ISIGN (1, IN-I) 

DJ = ISIGN (1,JN-J) 

DK - ISIGN (1,KN-K) 

I -IN 

J - JN 

K = KN 

GOTO 20 

C 

C We've been here before... 

C 

ELSE 

C 

C It seems to be outside the domain. We would have to extrapolate to 
C find it. 

C 

STATUS- 1 

IF (DIALOG) WRITE (6,*) 'EXTRAPOLATE' 

END IF 
END IF 
C 

50 CONTINUE 



RETURN 

END 

*** End of Subroutine Cell3 


★ ★****************************************************** ******** ******** 
SUBROUTINE CL0SE3 { ISRCH, I GRID, NSUBS, I, J, K, A, B, G, 

C X, Y, Z, I SUB, STATUS) 

a*********************************************************************** 

C 

C Find the cell containing the given (x,y,z) values and return its 
C (i,j,k) indices and subset number. We will use several different 
C strategies to try to do this efficiently, based on the value of 
C ISRCH : 

C 

C ISRCH* 0 - General search. 

C 1 - Search from previous point. 

C -n - Coming from grid n. 

C 

C STATUS=1 - (x, y , z) point could not be found. 

C 

common /grdxyz/ xyz (151, 40, 110, 3), idim, jdim, kdim, 

& iblank (151,40,110), blank, isubs (2) , 

& jsubs (2) , ksubs (2 ) 

common /save/ igsav, issav, isav, jsav, ksav, asav, bsav, gsav, 

& first, back, second, ifatal 

LOGICAL BLANK, first, second, back 
• INTEGER STATUS 
C 

STATUS* 0 
C 

C Search from previous point. 

C 

IF ( ISRCH. EQ.l) THEN 

IF (IGRID.EQ. IGSAV .AND. ISSAV. GE . 1 .AND. ISSAV . LE .NSUBS ) THEN 
I SUB * ISSAV 

I * ISAV 

J - JSAV 

K * KSAV 

A = ASAV 

B * BSAV 

G = GSAV 

is = isubs (1) 
ie = isubs (2) 
js = jsubs (1) 
je * jsubs (2) 
ks * ksubs (1) 
ke = ksubs (2) 

CALL PSRCH3 (is, ie, js, je, ks, ke, I, J, K, A, B, G, X, Y, Z, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 
ENDIF 
C 

C Try edges next . 

C 

CALL ESRCH3 ( I , J, K, A, B, G, X, Y, Z, ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 
C 

C Try faces next . 

C 

CALL FSRCH3 ( I , J, K, A, B, G, X, Y, Z, ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 
C 

C Try hole boundaries next. 

C 
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CALL HSRCH3 ( I, J, K, A, B, G, X, Y, Z, ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 

Coming from grid n... 

ELSE IF (ISRCH.LT.O) THEN 

Try edges first . 

CALL ESRCH3 ( I , J, K, A, B, G, X, Y, Z , ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 

If that fails, try closest point with I BLANK ■ -ISRCH. 

CALL NSRCH3 (-ISRCH, I, J, K, A, B, G,X, Y, Z, ISUB, ISTAT) 
IF (ISTAT. EQ.O) GOTO 10 

General search. 

ELSE 

Try edges . 

CALL ESRCH3 ( I, J, K, A, B, G, X, Y, Z, ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 

Try faces next . 

CALL FSRCH3 ( I, J, K, A, B, G, X, Y, Z, ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 

Try hole boundaries next . 

CALL HSRCH3 (I,J,K,A, B,G,X,Y,Z, ISUB, ISTAT) 

IF (ISTAT. EQ.O) GOTO 10 
END IF 

Nothing worked. Sorry. 

STATUS= 1 
IGSAV - 0 
C 

10 CONTINUE 

IGSAV = I GRID 
ISSAV = ISUB 

ISAV = I 

JSAV = J 

KSAV = K 

ASAV = A 

BSAV - B 

GSAV - G 

RETURN 
END 

*** End of Subroutine Close3 


************************************************************************ 

SUBROUTINE DMIN3 ( IS, IE, JS, JE, KS, KE, X, Y, Z,MPTS, NPTS, IPTS,D2MIN) 
************************************************************************ 

C 

C Add to a list of minimum distance points <i,j,k) to (x,y,z). Note 

C that the calling program is responsible for maintaining D2MIN if the 

C IPTS list is to be added to. Don't overflow IPTS, but don't worry 
C about a message or anything. 
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common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151, 40, 110) , blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

DIMENSION IPTS (3,MPTS) 

LOGICAL BLANK 

No previous points. Establish a minimum distance. 

IF (NPTS.EQ.O) D2MIN= 1.E35 
DO 10 K- KS,KE 
DO 10 J- JS,JE 
DO 10 I- IS, IE 

IF (.NOT. BLANK .OR. IBLANK (I, J,K) .NE . 0) THEN 
D2= (X— XYZ (I, J,K, 1) ) **2+ (Y-XYZ (I, J,K,2) ) **2 
C + (Z-XYZ (I,J,K, 3))**2 

Got a closer point. This replaces all previous points. 

IF (D2 .LT.D2MIN) THEN 
D2MIN= D2 
NPTS = 1 
IPTS ( 1 , NPTS ) = I 
IPTS (2, NPTS )= J 
IPTS <3, NPTS) “ K 

Tie for closest point. Add this to the list. 

ELSE IF (D2 . EQ.D2MIN) THEN 
IF (NPTS.LT.MPTS) THEN 
NPTS= NPTS+1 
IPTS (1, NPTS) = I 
IPTS (2,NPTS)= J 
IPTS (3, NPTS) = K 
END IF 
END IF 
ENDIF 

10 CONTINUE 
RETURN 
END 

* End of Subroutine DMin3 


***************************** ****************** ************************ 

SUBROUTINE ESRCH3 ( I , J, K, A, B, G, X, Y, Z, I SUB, STATUS) 

************* ************************************************ ********** 

Search from closest point on any edge. 

STATUS=1 - Couldn't find the point. 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151, 40,110), blank, isubs (2) , 

& jsubs (2) , ksubs (2) 

LOGICAL BLANK 
INTEGER STATUS 

PARAMETER (MPTS=5) 

DIMENSION IPTS (3, MPTS) 

STATUS= 0 

Loop through the subsets . 



IS- ISUBS(l) 

IE- ISUBS (2) 

JS= JSUBS(l) 

JE= JSUBS (2) 

KS- KSUBS(l) 

KE— KSUBS (2) 

C 

C Find the closest point (s) on any edge. Be careful of degenerate 
C subsets . 

C 

NETS - 0 
D2MIN- 1.E35 
C 

CALL DMIN3 (IS, IE, JS, JS, KS,KS, X, Y, Z,MPTS, NPTS , IPTS, D2MIN) 

IF (JS.NE.JE) THEN 

CALL DMIN3 (IS, IE, JE, JE, KS, KS, X, Y, Z,MPTS, NPTS, IPTS,D2MIN) 

END IF 

IF (KS.NE.KE) THEN 

CALL DMIN3 (IS, IE, JS, JS, KE, KE, X, Y, Z,MPTS, NPTS, IPTS,D2MIN) 

IF (JS.NE.JE) THEN 

CALL DMIN3 (IS, IE, JE, JE, KE, KE, X, Y, Z,MPTS, NPTS, IPTS,D2MIN) 
END IF 
END IF 
C 

CALL DMIN3 (IS, IS, JS+1, JE-1, KS, KS, X, Y, Z,MPTS,NPTS, IPTS,D2MIN) 

IF (IS.NE.IE) THEN 

CALL DMIN3 ( IE , IE , JS+1 , JE- 1 , KS , KS , X, Y, Z , MPTS, NPTS , IPTS , D2MIN) 
END IF 

IF (KS.NE.KE) THEN 

CALL DMIN3 ( IS , IS , JS+1 , JE-1 , KE, KE, X, Y, Z, MPTS , NPTS , IPTS , D2MIN) 
IF (IS.NE.IE) THEN 

CALL DMIN3 (IE, IE, JS+1, JE-1,KE,KE, 

C X,Y, Z, MPTS, NPTS, IPTS, D2MIN) 

END IF 
END IF 
C 

CALL DMIN3 ( IS, IS, JS, JS, KS+1, KE-1, 

C X, Y, Z, MPTS, NPTS, IPTS, D2MIN) 

IF (IS.NE.IE) THEN 

CALL DMIN3(IE,IE, JS, JS, KS+1, KE-1, 

C X,Y, Z, MPTS, NPTS, IPTS, D2MIN) 

END IF 

IF (JS.NE.JE) THEN 

CALL DMIN3 (IS, IS, JE, JE, KS+1, KE-1, 

C X,Y, Z, MPTS, NPTS, IPTS, D2MIN) 

IF (IS.NE.IE) THEN 

CALL DMIN3 ( IE, IE, JE, JE, KS+1 , KE— 1, 

C X, Y,Z, MPTS, NPTS, IPTS, D2MIN) 

END IF 
END IF 
C 

C Search from each of these points . 

C 

DO 10 11= 1 , NPTS 
C I SUB = I I SUB 

isub = 1 

I = IPTS (1,11) 

J = IPTS (2,11) 

K = IPTS (3, II) 

A =0. 

B =0. 

G =0. 

CALL PSRCH3 (is, ie, js, je,ks,ke, I, J, K, A, B,G, X, Y, Z, ISTAT) 

IF ( ISTAT. EQ.O) GOTO 30 
10 CONTINUE 

C 
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Failed all searches . 

STATUS= 1 

30 CONTINUE 
RETURN 
END 

** End of Subroutine ESrch3 


*********************************************************************** 
SUBROUTINE FSRCH3 (I, J, K, A, B,G, X, Y, Z, ISUB, STATUS) 


Search from closest point on any face. 

STATUS-1 - Couldn't find the point. 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151, 40, 110) , blank, isubs (2) , 

& jsubs (2) ,ksubs (2) 

LOGICAL BLANK 
INTEGER STATUS 

PARAMETER (MPTS=5) 

DIMENSION IPTS (3,MPTS) 

STATUS” 0 

Loop through the subsets. 

IS= ISUBS (1) 

IE= ISUBS (2) 

JS= JSUBS (1) 

JE= JSUBS (2) 

KS= KSUBS(l) 

KE= KSUBS (2 ) 

Find the closest point (s) on any face. Be careful of degenerate 
subsets . 

NPTS = 0 
D2MIN” 1.E35 

CALL DMIN3 (IS, IE, JS, JE, KS, KS, 

X, Y, Z,MPTS, NPTS, IPTS,D2MIN) 

IF (KS.NE.KE) THEN 

CALL DMIN3 (IS, IE, JS, JE,KE,KE, 

X, Y, Z,MPTS , NPTS, IPTS,D2MIN) 

CALL DMIN3 (IS, IE, JS, JS , KS+1, KE-1, 

X, Y, Z, MPTS , NPTS, IPTS,D2MIN) 

IF (JS.NE.JE) THEN 

CALL DMIN3 (IS, IE, JE, JE, KS + 1, KE— 1, 

X,Y,Z, MPTS, NPTS, IPTS, D2MIN) 

END IF 

CALL DMIN3 (IS, IS, JS+1, JE-1, KS+1, KE-1, 

X, Y,Z, MPTS, NPTS, IPTS, D2MIN) 

IF (IS.NE.IE) THEN 

CALL DMIN3 (IE, IE, JS+1 , JE-1, KS+1, KE-1, 

X, Y, Z, MPTS, NPTS, IPTS, D2MIN) 

END IF 
END IF 

Search from each of these points. 


C 

C 

C 

C 

C 

C 
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c 

DO 10 II- 1,NPTS 
C I SUB - I I SUB 

isub — 1 

I = IPTS (1,11) 

J - IPTS (2 r II) 

K = IPTS (3 f II) 

A =0. 

B =0. 

G - 0. 

CALL PSRCH3 (IS, IE, JS, JE, KS, KE, 

C I, J, K, A, B, G, X, Y, Z, ISTAT) 

IF ( ISTAT. EQ.O) GOTO 30 
10 CONTINUE 

Failed all searches . 

STATUS- 1 

30 CONTINUE 
RETURN 
END 

** End of Subroutine FSrch3 


*********************************************************************** 

SUBROUTINE HSRCH3 ( I , J, K, A, B, G, X, Y, Z, ISUB, STATUS) 
*********************************************************************** 

Search from closest point on a hole boundary. 

STATUS-1 - Couldn't find the point. 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank (151, 40, 110) , blank, isubs (2) , 

& jsubs (2) , ksubs (2) 

LOGICAL BLANK 
INTEGER STATUS 

PARAMETER (MPTS-5) 

DIMENSION IPTS (3,MPTS) 

STATUS- 0 

Loop through the subsets . 

if (blank) then 

IS- ISUBS (1) 

IE- ISUBS (2) 

JS- JSUBS (1) 

JE- JSUBS (2) 

KS- KSUBS (1) 

KE- KSUBS (2) 

Search for the closest point (s) on the boundary of any hole. 

NPTS = 0 
D2MIN- 1.E35 
DO 10 KK- KS, KE 
DO 10 JJ- JS, JE 
DO 10 II- IS, IE 

IF (IBLANK(II, JJ,KK) .NE.O .AND. 

C ( (II.NE.l .AND. IBLANK(II-1, JJ,KK) .EQ.O) 

C .OR. (II. NE. IDIM .AND. IBLANK (I 1+1, JJ,KK) .EQ. 0) 



C .OR. (JJ.NE.l .AND. IBLANK (II, JJ-1,KK) .EQ . 0 ) 

C .OR. (JJ.NE.JDIM .AND. IBLANK (II, JJ+1,KK) .EQ . 0) 

C .OR. (KK.NE.l .AND. IBLANK ( II, JJ, KK-1) .EQ . 0 ) 

C .OR. (KK.NE.KDIM .AND. IBLANK (II, JJ f KK+1) .EQ . 0) ) ) THEN 

D2= (X— XYZ ( II , JJ, KK, 1) ) **2 

C +(Y-XYZ(II, JJ,KK,2) ) **2 

C + (Z-XYZ (II, JJ,KK, 3) ) **2 

IF (D2 .LT.D2MIN) THEN 
D2MIN - D2 

NPTS = 1 

IPTS(1,1)= II 
IPTS (2,1) = JJ 
IPTS (3, 1)- KK 
C 

C If several points have the same minimum distance, keep track of them 
C Use PSRCH3 to determine which is correct . 

C 

ELSE IF (D2.EQ.D2MIN) THEN 
IF (NPTS.LT.MPTS) THEN 
NPTS = NPTS+1 

IPTS (1,NPTS)= II 
IPTS (2, NPTS) = JJ 
IPTS (3, NPTS) = KK 
END IF 
END IF 
END IF 

10 CONTINUE 

C 

C Check each of these points with PSRCH3 . 

C 

DO 20 11= 1 , NPTS 
I SUB = 1 

I = IPTS (1,11) 

J = IPTS (2, II) 

K = IPTS (3, II) 

A = 0 . 

B =0. 

G — 0 . 

CALL PSRCH3 (IS, IE, JS, JE,KS,KE, 

C I,J,K,A,B,G,X,Y,Z, ISTAT) 

IF ( ISTAT. EQ.0) GOTO 30 
20 CONTINUE 

C 

ENDIF 

C 

C Failed from all closest hole boundary points . 

C 

STATUS= 1 
C 

30 CONTINUE 
RETURN 
END 

*** End of Subroutine HSrch3 


************************************************************************ 
SUBROUTINE INV3X3 (A, AINV, STATUS) 

★★A**************************** ★*★★**★**★*★★**★*★**** *★*★***★*★*★★*★**** 
C 

C Invert the 3x3 matrix A. If A is singular, do our best to find the 
C pseudo-inverse. 

C 

C STATUS=1 - A has one dependent column. 

C STATUS=2 - A has two dependent columns . 

C STATUS=3 - A is zero. 



c 


c 

c 

c 


c 


c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


DIMENSION A(3,3) ,AINV(3,3) 

INTEGER STATUS 

DIMENSION TMP (3,3) , WORK (3) , S (3) , E (3) ,U (3, 3) ,V(3,3) ,SIU(3, 3) 
STATUS- 0 

AINV (1, 1) — A(2,2) *A(3,3) - A (3, 2) *A(2, 3) 

AINV(1, 2) — A (3, 2) *A(1,3) - A <1, 2) *A<3, 3) 

AINV (1,3)— A(l,2) *A(2, 3) - A (2, 2) *A (1, 3) 

AINV (2,1)= A(3, 1) *A (2, 3) - A(2, 1) *A(3, 3) 

AINV (2,2)— A(l, 1) *A (3,3) - A{3, 1) *A(1,3) 

AINV (2, 3) — A(2, 1) *A(1, 3) - A(l, 1) *A(2, 3) 

AINV (3, 1) = A(2, 1) *A(3,2) - A(3, 1) *A(2,2) 

AINV (3,2) — A(3, 1) *A(1,2) - A(l, 1) *A(3,2) 

AINV (3,3) = A (1, 1 ) *A (2, 2 ) - A(2,1)*A(1,2) 

DET— A(l, 1) *AINV(1, 1) + A(2, 1) *AINV(1,2) + A(3, 1) *AINV(1, 3) 

Matrix is nonsingular. Finish up AINV. 

IF (DET.NE.O.) THEN 
DET- 1 . /DET 

AINV ( 1 , 1 ) — AINV (1,1) *DET 
AINV (2,1)= AINV (2,1) *DET 
AINV (3,1)= AINV (3,1) *DET 
AINV (1,2)= AINV (1,2) *DET 
AINV (2,2)= AINV (2,2) *DET 
AINV (3,2)= AINV (3,2) *DET 
AINV (1,3)— AINV (1,3) *DET 
AINV (2,3)= AINV (2,3) *DET 
AINV (3,3)= AINV (3,3) *DET 

Matrix is singular. Do a singular value decomposition to construct 
the pseudo-inverse. Use UNPACK routine SSVDC. 

ELSE 

CALL COPY (9, A, TMP) 

CALL SSVDC (TMP, 3,3,3,S,E,U,3,V, 3, WORK, 11, INFO) 

IF (S(l).EQ.O.) THEN 
STATUS- 3 
CALL ZERO (9, AINV) 

GOTO 10 
ENDIF 

-1 

Compute VS U. 

S(l)= l./S(l) 

IF (S(3)*S(1) .LT.l.E-5) THEN 
STATUS- 1 
S (3) =0. 

ELSE 

S (3) = l./S(3) 

ENDIF 

IF (S(2)*S(1) .LT.l.E-5) THEN 
STATUS- 2 
S (2) =0. 

ELSE 

S (2) - l./S(2) 

ENDIF 

Start out assuming S is a diagonal matrix. 
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SIU<1,1) - S(1)*U(1,1) 

SIU (2,1) - S (2) *U (2, 1) 

SIU (3,1) = S (3) *U (3, 1) 

SIU(1,2) - S (1) *U ( 1, 2) 

SIU (2,2) - S (2) *U (2,2) 

SIU (3, 2) - S(3)*U(3,2) 

SIU (1, 3) = S(1)*U(1,3) 

SIU (2, 3) - S (2) *U (2, 3) 

SIU (3, 3) - S (3) *U (3, 3) 

S is upper bidiagonal, with E as the super diagonal. 

IF (INF0.GE.1) THEN 

SIU(1,1) = SIU(1,1) - (E(l) *S (1) *S (2) ) *U(2,1) 

SIU (1, 2) = SIU (1,2) - (E (1) *S (1) *S (2) ) *U(2,2) 

SIU (1,3) = SIU(1,3) - (E (1) *S (1) *S (2) ) *U(2, 3) 

END IF 

IF (INFO.GE.2) THEN 

SIU(1,1) = SIU(1,1) + (E(l) *E(2) *S(1) *S (2) **2*S (3) ) *U (3, 1) 
SIU (1,2) = SIU(1,2) + (E(1)*E(2)*S(1)*S(2)**2*S(3))*U(3,2) 
SIU (1, 3) = SIU(1,3) + (E(1)*E(2)*S(1)*S(2)**2*S(3))*U(3,3) 

SIU (2, 1) = SIU (2, 1) - (E(2)*S(2)*S(3))*U(3,1) 

SIU (2, 2) - SIU (2,2) - (E (2) *S (2) *S (3) ) *U(3,2) 

SIU (2, 3) - SIU(2, 3) - (E(2) *S (2) *S (3) ) *U(3,3) 

ENDIF 

+ -1 

Finish up A = V S U. 

AINV (1,1)= V(1,1)*SIU(1,1) + V (1, 2) *SIU (2,1) + V(l, 3) *SIU (3, 1) 
AINV (2, 1) = V(2, 1) *SIU (1,1) + V(2, 2) *SIU (2,1) + V (2, 3) *SIU <3, 1) 
AINV(3,1)= V(3, 1) *SIU (1, 1) + V(3,2)*SIU(2,1) + V(3, 3) *SIU (3, 1) 

AINV(1,2)= V(1,1)*SIU(1,2) + V(1,2)*SIU<2,2) + V<1, 3) *SIU (3, 2) 
AINV (2,2)= V(2, 1) *SIU (1,2) + V(2,2) *SIU(2,2) + V (2, 3) *SIU (3, 2) 
AINV (3,2)= V(3, 1) *SIU (1,2) + V (3,2) *SIU(2,2) + V{3, 3) *SIU (3, 2 ) 

AINV (1,3)= V(1,1)*SIU(1,3) + V (1, 2) *SIU (2,3) + V(l, 3) *SIU (3, 3) 
AINV (2,3)= V(2, 1) *SIU (1,3) + V (2, 2) *SIU (2, 3) + V (2, 3) *SIU (3, 3) 
AINV (3,3)= V (3, 1) *SIU ( 1 , 3 ) + V(3,2)*SIU(2,3) + V(3, 3) *SIU (3, 3) 
ENDIF 

10 CONTINUE 
RETURN 
END 

End of Subroutine Inv3X3 


SUBROUTINE NSRCH3 (N, I, J, K, A, B, G, X, Y, Z, ISUB, STATUS) 


Search from closest point coming from grid N. 

STATUS=1 - Couldn't find the point. 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 
i iblank (151, 40, 110) , blank, isubs (2) , 

& jsubs (2) , ksubs (2) 

LOGICAL BLANK 



non oooo ooo ooo 


INTEGER STATUS 
C 

PARAMETER (MPTS=5) 

DIMENSION IPTS (3, MPTS) 

C 

STATUS= 0 

Loop through the subsets . 

IF (BLANK) THEN 

IS- ISUBS(l) 

IE- ISUBS (2) 

JS- JSUBS(l) 

JE- JSUBS (2) 

KS- KSUBS(l) 

KE- KSUBS (2) 

Search for the closest point (s) on the boundary of any hole. 

NPTS - 0 
D2MIN- 1.E35 
DO 10 KK= KS, KE 
DO 10 JJ= JS, JE 
DO 10 11= IS, IE 

IF (IBLANKdl, JJ,KK) .EQ.-N) THEN 
D2— (X-XYZ (II, JJ , KK, 1) ) **2 
C +(Y-XYZ(II, JJ,KK,2) ) **2 

C + (Z-XYZ (II, JJ, KK, 3) ) **2 

IF (D2.LT.D2MIN) THEN 
D2MIN - D2 

NPTS = 1 

IPTS (1,1)— II 
IPTS (2,1)= JJ 
IPTS (3,1)= KK 

If several points have the same minimum distance, keep track of them. 
Use PSRCH3 to determine which is correct . 

ELSE IF (D2.EQ.D2MIN) THEN 
IF (NPTS. LT. MPTS) THEN 
NPTS = NPTS+1 

IPTS (1, NPTS) = II 
IPTS (2, NPTS) = JJ 
IPTS (3, NPTS) = KK 
ENDIF 
ENDIF 
ENDIF 

10 CONTINUE 

Check each of these points with PSRCH3 . 

DO 20 11= 1, NPTS 
I SUB = 1 

I = IPTS (1, II) 

J = IPTS (2, II) 

K = IPTS (3, II) 

A =0. 

B =0. 

G =0. 

CALL PSRCH3 (IS, IE, JS, JE, KS, KE, 

C I,J,K,A, B,G,X,Y,Z, ISTAT) 

IF ( ISTAT. EQ.0) GOTO 30 
20 CONTINUE 

C 


C 


ENDIF 



C Failed from all closest hole boundary points. 
C 

STATUS* 1 
C 

30 CONTINUE 
RETURN 
END 

*** End of Subroutine NSrch3 


**************************************************************** ******** 

SUBROUTINE PSRCH3 (IS, IE, JS, JE, KS, KE, I, J, K, A, B, G, X, Y, Z, STATUS ) 
★★a********************************************************** *********** 

C 

C Search for a point using CELL3 • 

C 

C STATUS=1 - Couldn't find the point. 

C 

common /grdxyz/ xyz (151, 40, 110, 3) , idim, jdim, kdim, 

& iblank(151, 40, 110) , blank, isubs (2) , 

& jsubs (2) , ksu bs (2) 

LOGICAL BLANK 
INTEGER STATUS 
C 

DIMENSION AMAT (3,3) , BMAT (3,3) 

C 

STATUS* 0 
C 

CALL CELL3 (IS , IE, JS , JE, KS, KE, 

C I,J,K,A,B,G,X,Y, 2, AMAT, BMAT, ISTAT) 

IF ( ISTAT. NE.0) THEN 
STATUS* 1 
GOTO 10 
END IF 
C 

10 CONTINUE 
RETURN 
END 

*** End of Subroutine PSrch3 


************************************************************* *********** 
SUBROUTINE COPY (LEN, FROM, TO) 

********************* ****************************** ********************* 
C 

c Copy an array from FROM to TO. The user is responsible for making 
C sure that FROM doesn't overwrite itself if FROM and TO overlap. 

C 

DIMENSION FROM ( * ) , TO ( * ) 

C 

DO 10 1= 1, LEN 
TO ( I) = FROM (I) 

10 CONTINUE 
RETURN 
END 

*** End of Subroutine Copy 


************************************************************************ 
subroutine ssvdc (x, ldx, n,p,s,e,u, ldu, v, ldv, work, job, info) 



********** ************************************************************** 
integer ldx,n,p, ldu, ldv, job, info 

real x (ldx, 1) , s (1) , e (1) ,u (ldu, 1) , v(ldv, 1) , work (1) 
c 

Caveat receptor. (Jack) dongarra@anl-mcs, (Eric Grosse) research !ehg 
Compliments of netlib Fri Dec 12 16:40:31 CST 1986 
c 

c ssvdc is a subroutine to reduce a real nxp matrix x by 

orthogonal transformations u and v to diagonal form, the 
diagonal elements s(i) are the singular values of x. the 
columns of u are the corresponding left singular vectors, 
and the columns of v the right singular vectors. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


on entry 


ldx 


ldu 


ldv 


work 


job 


on return 


real ( ldx, p) , where ldx.ge.n. 
x contains the matrix whose singular value 
decomposition is to be computed, x is 
destroyed by ssvdc. 

integer . 

ldx is the leading dimension of the array x. 


integer . 
n is the 

integer, 
p is the 


number of rows of the matrix x. 


number of columns of the matrix x. 


integer . 

ldu is the leading dimension of the array u. 
(see below) . 

integer. 

ldv is the leading dimension of the array v. 
(see below) . 


real (n) 
work is 


a scratch array. 


integer. 

job controls the computation of the singular 
vectors. it has the decimal expansion ab 
with the following meaning 


a .eq. 0 
a .eq. 1 

a . ge . 2 

b. eq. 0 
b.eq. 1 


do not compute the left singular 
vectors . 

return the n left singular vectors 
in u. 

return the first min(n,p) singular 
vectors in u. 

do not compute the right singular 
vectors. 

return the right singular vectors 
in v. 


real (mm) , where mm=min (n+1, p) . 
the first min(n,p) entries of s contain the 
singular values of x arranged in descending 
order of magnitude. 

real (p) . 

e ordinarily contains zeros, however see the 
discussion of info for exceptions. 



u 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


real (ldu,k) , where ldu.ge.n. if joba.eq.l then 
k.eq.n, if joba.ge.2 then 
k.eq.min (n, p) . 

u contains the matrix of left singular vectors, 
u is not referenced if joba.eq.O. if n.le.p 
or if joba.eq.2, then u may be identified with x 
in the subroutine call. 

real(ldv,p), where ldv.ge.p. 

v contains the matrix of right singular vectors . 
v is not referenced if job.eq.O. if p.le.n, 
then v may be identified with x in the 
subroutine call. 


info integer. 

the singular values (and their corresponding 
singular vectors) s (info+1) , s (info+2) , . . . , s (m) 
are correct (here m=min{n,p)). thus if 
info. eq. 0, all the singular values and their 
vectors are correct. in any event, the matrix 
b = trans(u)*x*v is the bidiagonal matrix 
with the elements of s on its diagonal and the 
elements of e on its super-diagonal (trans (u) 
is the transpose of u) . thus the singular 
values of x and b are the same. 

Unpack, this version dated 03/19/79 . 

correction to shift calculation made 2/85. 

g.w. Stewart, university of maryland, argonne national lab. 

***** uses the following functions and subprograms . 


external srot 

bias saxpy, sdot, sscal, sswap, snrm2, srotg 
fortran abs, amaxl,maxO,minO,mod, sqrt 


internal variables 


integer i, iter, j , jobu, k, kase, kk, 1, 11, 11s, 1ml, lpl, Is, lu,m,maxit, 
mm, mml , mpl , net , nctpl , ncu, nrt , nrtpl 
real sdot,t,r 

real b, c, cs, el, emml , f , g, snrm2, scale, shift, si, sm, sn, smml, tl, test , 
ztest 

logical wantu,wantv 


set the maximum number of iterations, 
maxit = 30 

determine what is to be computed. 

wantu = .false, 
wantv = .false, 
jobu = mod( job, 100) /10 
ncu * n 

if (jobu .gt. 1) ncu = min0(n,p) 
if (jobu ,ne. 0) wantu = .true, 
if (mod (job, 10) .ne. 0) wantv = .true. 

reduce x to bidiagonal form, storing the diagonal elements 
in s and the super-diagonal elements in e. 

info = 0 

net = min0(n-l,p) 



c 

c 

c 

c 


10 

20 


c 

c 

c 


30 


c 

c 

c 

c 


40 

50 


c 

c 

c 

c 


60 

70 


c 

c 

c 

c 


80 


c 

c 

c 


90 

100 


nrt * maxO (0,min0 (p-2,n) ) 
lu * maxO (net, nrt) 
if (lu .It. 1) go to 170 
do 160 1 * 1, lu 
Ipl - 1 + 1 

if (1 .gt. net) go to 20 

compute the transformation for the 1-th column and 
place the 1-th diagonal in s(l). 

s(l) - snrm2 (n-1+1, x (1,1) , 1) 
if (s(l) .eq. O.OeO) go to 10 

if (x (1, 1) .ne. O.OeO) s (1) « sign(s (1) ,x(l, 1) ) 
call sscal (n-1+1 , l.OeO/s (1) ,1) 
x (1, 1) = l.OeO + x(U) 
continue 
3(1) = -s(l) 
continue 

if (p .It. lpl) go to 50 
do 40 j 2 lpl# p 

if (1 .gt. net) go to 30 
if (s(l) .eq. O.OeO) go to 30 

apply the transformation. 

t = -sdot (n-1+1, x(l,l) ,l,x(l, j),l)/x(l,l) 
call saxpy (n-1+1, t, x (1,1) , l,x(l, j) , 1) 
continue 

place the 1-th row of x into e for the 
subsequent calculation of the row transformation. 

e(j) = x(l,j) 
continue 
continue 

if (.not.wantu .or. 1 .gt. net) go to 70 

place the transformation in u for subsequent back 
multiplication . 

do 60 i = 1, n 

u (i r 1) = x (i, 1) 
continue 
continue 

if (1 .gt. nrt) go to 150 

compute the 1-th row transformation and place the 
1-th super-diagonal in e(l). 

e(l) = snrm2 (p-l,e (lpl) , 1) 
if (e(l) .eq. O.OeO) go to 80 

if (e (lpl) .ne. O.OeO) e(l) - sign <e (1) ,e dpi) ) 
call sscal (p-1, 1 . OeO/e (1) ,e(lpl) ,1) 
e(lpl) = l.OeO + e(lpl) 
continue 
e (1) « -e(l) 

if (lpl .gt. n .or. e(l) .eq. O.OeO) go to 120 

apply the transformation. 

do 90 i = lpl/ n 
work(i) -= O.OeO 
continue 

do 100 j = lpl, p 

call saxpy (n-1, e ( j ) , x ( lpl , j ) , 1 , work ( lpl ) , 1 ) 
continue 



110 

120 


do 110 j = lpl, p 

call saxpy (n-l,-e ( j) /e (lpl) , work (lpl) , l,x(lpl, j) , 1) 
continue 
continue 

if (.not.wantv) go to 140 
c 

c place the transformation in v for subsequent 

c back multiplication, 

c 

do 130 i * lpl, p 
v(i,l) - e (i) 

130 continue 

140 continue 

150 continue 
160 continue 
170 continue 


c 

c set up the final bidiagonal matrix or order m. 

c 

m = minO (p, n+1) 
nctpl = net + 1 
nrtpl = nrt + 1 

if (net .It. p) s (nctpl) * x (nctpl, nctpl) 
if <n .It. m) s (m) = O.OeO 
if (nrtpl .It. m) e (nrtpl) = x (nrtpl, m) 
e(m) = O.OeO 


c 

c 

c 


if 

if 


180 

190 

200 


210 

220 


230 

240 

250 


260 

270 

280 

290 


required, generate u. 

(.not.wantu) go to 300 
if (ncu .It. nctpl) go to 200 
do 190 j * nctpl, ncu 
do 180 i * 1, n 
u (i, j) = 0 . OeO 
continue 
u(j,j) = l.OeO 
continue 
continue 

if (net .It. 1) go to 290 
do 280 11 = 1, net 
1 = net - 11 + 1 
if (s (1) .eq. O.OeO) go to 250 
lpl = 1 + 1 

if (ncu .It. lpl) go to 220 
do 210 j = lpl, ncu 

t = -sdot (n-l+l,u (1, 1) ,l,u(l, j) , 1) /u(l,l) 
call saxpy (n-l+l,t,u (1, 1) ,l,u(l, j) ,1) 
continue 
continue 

call sscal (n-1+1, -1 .0e0,u(l,l),l) 
u(l,l) = l.OeO + u(l,l) 
lml =1-1 

if (1ml .It. 1) go to 240 
do 230 i = 1, lml 
u (i, 1) = 0 . OeO 
continue 
continue 
go to 270 
continue 

do 260 i * 1, n 
u(i,l) = O.OeO 
continue 
u(l,l) = l.OeO 
continue 
continue 
continue 



300 continue 


c 

c if it is required, generate v. 

c 

if (.not.wantv) go to 350 
do 340 11 = 1, p 

1 * p - 11 + 1 

lpl - 1 + 1 

if (1 .gt. nrt ) go to 320 
if (e (1) .eq. O.OeO) go to 320 
do 310 j =■ lpl, p 

t * -sdot (p-1, v (lpl, 1) , l,v(lpl, j) ,1) /v(lpl,l) 
call saxpy (p-1, t , v (lpl, 1) , 1, v (lpl, j) ,1) 

310 continue 

320 continue 

do 330 i * 1, p 
v(i,l) « O.OeO 
330 continue 

v ( 1 , 1 > = l.OeO 
340 continue 
350 continue 
c 

c main iteration loop for the singular values, 

c 

mm — m 
iter * 0 


360 continue 
c 

c quit if all the singular values have been found, 

c 

c . . .exit 

if (m .eq. 0) go to 620 
c 

c if too many iterations have been performed, set 

c flag and return, 

c 

if (iter .It. maxit) go to 370 
info = m 

c exit 

go to 620 
370 continue 
c 

c this section of the program inspects for 

c negligible elements in the s and e arrays, on 

c completion the variables kase and 1 are set as follows, 

c 

c kase =1 if s (m) and e(l-l) are negligible and l.lt.m 

c kase =2 if s(l) is negligible and l.lt.m 

c kase * 3 if e(l-l) is negligible, l.lt.m, and 

c s(l) , s (m) are not negligible (qr step), 

c kase * 4 if e(m-l) is negligible (convergence), 

c 

do 390 11 = 1, m 
1 = m - 11 
c . . .exit 


if (1 .eq. 0) go to 400 
test = abs(s(l)) + abs(s(l+l)) 
ztest = test + abs(e(l)) 
if (ztest ,ne. test) go to 380 
e (1) = O.OeO 

c exit 

go to 400 
continue 
continue 
continue 

if (1 .ne. m - 1) go to 410 


380 

390 

400 



410 


c 


c 


420 

430 

440 


450 


460 


470 

480 


c 

c 

c 

c 

c 

c 

490 


500 

510 


c 

c 

c 


kase * 4 
go to 480 
continue 

lpl *1+1 
mpl * m + 1 
do 430 11s * lpl, mpl 
Is * m - 11s + lpl 
. . .exit 

if (Is .eq. 1) go to 440 
test * O.OeO 

if (Is .ne. m) test = test + abs(e(ls)) 
if (Is .ne. 1+1) test * test + abs(e(ls-l)) 
ztest = test + abs(s(ls)) 
if (ztest .ne. test) go to 420 
s(ls) = O.OeO 

exit 

go to 440 
continue 
continue 
continue 

if (Is .ne. 1) go to 450 
kase * 3 
go to 470 
continue 

if (Is .ne. m) go to 460 
kase = 1 
go to 470 
continue 
kase * 2 

I = Is 
continue 

continue 
1*1 + 1 

perform the task indicated by kase. 

goto (490,520,540,570), kase 

deflate negligible s (m) . 

continue 

mml = m - 1 
f * e(m-l) 
e(m-l) = O.OeO 
do 510 kk = 1, mml 
k = mml - kk + 1 
tl = s (k) 

call srotg (tl, f , cs, sn) 
s(k) = tl 

if (k .eq. 1) go to 500 
f = -sn*e (k-1) 
e(k-l) = cs*e (k-1) 
continue 

if (wantv) call srot (p, v (1, k) , 1, v (l,m) , 1, cs, sn) 
continue 
go to 610 

split at negligible s(l). 

continue 

f = e (1-1) 
e (1-1) * O.OeO 
do 530 k * 1, m 

I I = s (k) 

call srotg(tl,f,cs,sn) 
s (k) = tl 


520 



530 


f * -sn*e(k) 
e(k) - cs*e(k) 

if (wantu) call srot (n, u (l,k) , 1, u (1, 1-1) , 1, cs, sn) 
continue 
go to 610 

perform one qr step, 
continue 

calculate the shift. 

scale - amaxl (abs (s (m) ) , abs (s (m-1) ) , abs (e (m-1) ) , abs (s (1) ) , 
abs(e(l))) 
sm = s (m) /scale 
smml = s (m-1) /scale 
eittml = e (m-1) /scale 
si » s (1) /scale 
el * e (1) /scale 

b * ((smml + sm)*(smml - sm) + emml**2) /2 . OeO 
c = (sm*emml)**2 
shift = O.OeO 

if (b .eq. O.OeO .and. c .eq. O.OeO) go to 550 
shift * sqrt (b**2+c) 
if (b .It. O.OeO) shift * -shift 
shift = c/(b + shift) 
continue 

f = (si + sm)*(sl “ sm) + shift 
g = sl*el 

chase zeros . 

rami = m - 1 
do 560 k * 1, mml 

call srotg(f ,g,cs,sn) 
if (k .ne. 1) e(k-l) = f 
f - cs*s(k) + sn*e(k) 
e(k) = cs*e(k) - sn*s(k) 
g = sn*s (k+1) 
s (k+1) = cs*s(k+l) 

if (wantv) call srot (p, v (1, k) , 1, v (1, k+1) , 1, cs, sn) 
call srotg (f , g, cs, sn) 
s(k) * f 

f - cs*e (k) + sn*s (k+1) 
s(k+l) = -sn*e(k) + cs*s(k+l) 
g = sn*e(k+l) 
e (k+1) = cs*e (k+1) 
if (wantu .and. k .It. n) 

call srot (n,u(l,k),l,u(l,k+l),l,cs,sn) 
continue 
e (m-1) = f 
iter * iter + 1 
go to 610 

convergence . 

continue 

make the singular value positive. 

if (s (1) .ge. O.OeO) go to 580 
s(l) = -s (1) 

if (wantv) call sscal (p, -1 . OeO, v (1, 1) , 1) 
continue 


c 

c 


580 


order the singular value. 



c 


590 if (1 .eq. mm) go to 600 

c • exit 

’if (s(l) .ge. 3(1+1)) go to 600 
t - s(l) 
s(l) * s(l+l) 

3(1+1) * t 

if (wantv .and. 1 .It. p) 

* call sswap (p, v (1, 1) , 1, v (1, 1+1) , 1) 
if (wantu .and. 1 .It. n) 

* call sswap (n, u (1, 1) , l,u (1, 1+1) , 1) 
1-1 + 1 



go to 590 

600 

continue 


iter = 0 


m = m - 1 

610 

continue 


go to 360 
620 continue 
return 
end 

*** End of Subroutine SSVDC 


************************************************************************ 
subroutine sswap (n, sx, incx, sy, incy) 

★a********************************************************************** 

c 

c interchanges two vectors . 

c uses unrolled loops for increments equal to 1. 

c jack dongarra, Unpack, 3/11/78. 

c 

real sx ( 1 ) , sy ( 1 ) , stemp 
integer i, incx, incy, ix, iy,m,mpl, n 
c 

if (n . le . 0 ) return 

if (incx.eq. 1 .and. incy .eq. 1) go to 20 
c 

c code for unequal increments or equal increments not equal 

c to 1 

c 

ix - 1 
iy = 1 

if (incx . It . 0) ix - (-n+1) *incx + 1 
if (incy . It . 0) iy - (-n+1) *incy + 1 
do 10 i = l,n 
stemp * sx(ix) 
sx(ix) = sy(iy) 

sy(iy) = stemp 

ix = ix + incx 
iy = iy + incy 
10 continue 
return 
c 

c code for both increments equal to 1 

c 

c 

c clean-up loop 

c 

20 m = mod (n, 3) 

if ( m .eq. 0 ) go to 40 
do 30 i = l,m 
stemp = sx(i) 
sx(i) = sy (i) 



30 


sy{i) * steinp 
continue 

if ( n . It . 3 ) return 
40 mpl = m + 1 

do 50 i * mpl, n, 3 
steinp * sx(i) 
sx(i) - sy (i) 
sy(i) = stemp 
stomp = sx(i + 1) 
sx(i + 1) « sy(i + 1) 
sy(i + 1) - stemp 
steinp * sx(i + 2) 
sx(i + 2) - sy(i + 2) 
sy(i + 2) - stemp 
50 continue 
return 
end 

*** End of Subroutine SSwap 


************************************************************************ 
subroutine saxpy (n, sa, sx f incx, sy, incy) 
************************************************************************ 
c 

c constant times a vector plus a vector, 

c uses unrolled loop for increments equal to one. 
c jack dongarra, Unpack, 3/11/78. 
c 

real sx<l) , sy (1) , sa 
integer i, incx, incy, ix, iy,m,mpl, n 
c 

if (n . le . 0 ) return 
if (sa .eq. 0.0) return 
if ( incx. eq.l. and. incy. eq.l) go to 20 
c 

c code for unequal increments or equal increments 

c not equal to 1 

c 

ix = 1 

iy = 1 

if (incx.lt . 0) ix - (-n+l)*incx + 1 
if (incy.lt .0) iy = (-n+l)*incy + 1 
do 10 i = l,n 

sy(iy) = sy (iy) + sa*sx(ix) 
ix * ix + incx 

iy = iy t incy 

10 continue 
return 
c 

c code for both increments equal to 1 

c 

c 

c clean-up loop 

c 

20 m = mod (n, 4) 

if( m .eq. 0 ) go to 40 
do 30 i * l,m 

sy(i) = sy(i) + sa*sx(i) 

30 continue 

if ( n .It. 4 ) return 
40 mpl = m + 1 

do 50 i = mpl,n, 4 

sy(i) = sy (i) + sa*sx(i) 

sy (i + 1) = sy(i + 1) + sa*sx{i + 1) 
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sy(i + 2) « sy(i + 2) + sa*sx(i + 2) 
sy(i + 3) * sy(i + 3) + sa*sx(i + 3) 
50 continue 
return 
end 

*** End of Subroutine SaxPy 


******************** ************************************ ******* ********* 

real function sdot (n, sx, incx, sy, incy) 

*********************** ********************************** ****** ********* 

c 

forms the dot product of two vectors . 
uses unrolled loops for increments equal to one. 
c jack dongarra. Unpack, 3/11/78. 

c 

real sx ( 1 ) , sy ( 1 ) , stemp 
integer i, incx, incy/ ix, iy^/inpl, n 
c 

stemp = O.OeO 
sdot - O.OeO 
if (n . le . 0 ) return 

if (incx .eq. 1 . and. incy .eq. 1) go to 20 
c 

c code for unequal increments or equal increments 

c not equal to 1 

c 

ix * 1 
iy = 1 

if (incx. It . 0) ix = (-n+l)*incx + 1 
if (incy . It . 0) iy = (-n+l)*incy + 1 
do 10 i = l,n 

stemp = stemp + sx (ix) *sy (iy) 
ix = ix + incx 

iy = iy + incy 

10 continue 

sdot = stemp 
return 
c 

c code for both increments equal to 1 

c 

c 

c clean-up loop 

c 

20 m = mod (n, 5) 

if ( m .eq. 0 ) go to 40 
do 30 i = 1 , m 

stemp = stemp + sx(i)*sy(i) 

30 continue 

if ( n . It . 5 ) go to 60 
40 mpl = m + 1 

do 50 i = mpl,n f 5 

stemp = stemp + sx(i)*sy(i) + sx(i + l)*sy(i + 1) + 

* sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 
50 continue 
60 sdot = stemp 
return 
end 

*** End of Function SDot 
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real function snrm2 ( n, sx, incx) 

★ ★**★**★★*★★*★*****★*★*****★**★*★★★★**★★★*★★★*★*★★***•**★*★*★★**★*★★****★ 
integer next 

real sx(l) f cutlo, cuthi, hitest, sum, xmax, zero, one 
data zero, one /O.OeO, l.OeO/ 
c 

c euclidean norm of the n-vector stored in sx() with storage 

c increment incx . 

c if n . le . 0 return with result * 0 . 

c if n ,ge. 1 then incx must be .ge. 1 
c 

c c.l.lawson, 1978 jan 08 

c 

c four phase method using two built-in constants that are 

c hopefully applicable to all machines. 

c cut lo * maximum of sqrt(u/eps) over all known machines, 

c cuthi * minimum of sqrt (v) over all known machines, 

c where 

c eps * smallest no. such that eps + 1. ,gt. 1. 

c u * smallest positive no. (underflow limit) 

c v * largest no. (overflow limit) 

c 

c brief outline of algorithm. . 

c 

c phase 1 scans zero components . 

c move to phase 2 when a component is nonzero and .le. cutlo 

c move to phase 3 when a component is .gt. cutlo 

c move to phase 4 when a component is .ge. cuthi/m 

c where m = n for x() real and m = 2*n for complex, 

c 

c values for cutlo and cuthi.. 

c from the environmental parameters listed in the imsl converter 
c document the limiting values are as follows.. 

c cutlo, s.p. u/eps - 2** (-102) for honeywell. close seconds are 
c univac and dec at 2** (-103) 

c thus cutlo - 2** (-51) * 4 . 44089e-16 

c cuthi, s.p. v - 2**127 for univac, honeywell, and dec. 
c thus cuthi = 2** (63.5) * 1.30438el9 

c cutlo, d.p. u/eps = 2** (-67) for honeywell and dec. 
c thus cutlo = 2** (-33.5) = 8.23181d-ll 

c cuthi, d.p. same as s.p. cuthi ■ 1.30438dl9 
c data cutlo, cuthi / 8.232d-ll, 1.304dl9 / 

c data cutlo, cuthi / 4.441e-16, 1.304el9 / 

data cutlo, cuthi / 4.441e-16, 1.304el9 / 

c 

if (n .gt. 0) go to 10 
snrm2 * zero 
go to 300 
c 

10 assign 30 to next 
sum * zero 
nn « n * incx 

c begin main loop 

i = 1 

20 go to next, (30, 50, 70, 110) 

30 if ( abs(sx(i)) .gt. cutlo) go to 85 
assign 50 to next 
xmax = zero 
c 

c phase 1. sum is zero 

c 

50 if ( sx(i) .eq. zero) go to 200 

if ( abs(sx(i)) .gt. cutlo) go to 85 

prepare for phase 2. 

assign 70 to next 
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go to 105 

prepare for phase 4. 

100 i - j 

assign 110 to next 

sum - (sum / sx(i)) / sx(i) 

105 xmax * abs(sx(i)) 
go to 115 

phase 2. sum is small. 

scale to avoid destructive underflow. 

70 if ( abs(sx(i)) .gt. cutlo ) go to 75 

common code for phases 2 and 4 . 

in phase 4 sum is large, scale to avoid overflow. 

110 if ( abs(sx(i)) . le . xmax ) go to 115 

sum = one + sum * (xmax / sx(i))**2 
xmax * abs (sx (i) ) 
go to 200 

115 sum * sum + (sx{i) /xmax) **2 
go to 200 


c prepare for phase 3. 

c 

75 sum = (sum * xmax) * xmax 
c 
c 

c for real or d.p. set hitest = cuthi/n 

c for complex set hitest = cuthi/(2*n) 

c 

85 hitest = cuthi/float ( n ) 
c 

c phase 3. sum is mid-range, no scaling, 

c 

do 95 j =i,nn,incx 

if (abs (sx ( j ) ) ,ge. hitest) go to 100 
95 sum = sum + sx(j)**2 
snrm2 = sqrt ( sum ) 
go to 300 
c 

200 continue 

i = i + incx 

if ( i .le. nn ) go to 20 
c 

c end of main loop, 

c 

c compute square root and adjust for scaling, 

c 

snrm2 = xmax * sqrt (sum) 

300 continue 
return 
end 

*** End of Function SNRM2 


subroutine srot (n, sx, incx, sy, incy, c, s) 
c 



c applies a plane rotation, 

c jack dongarra, linpack, 3/11/78. 

c 

real sx (1) , sy (1) , stemp, c, s 
integer i, incx, incy, ix, iy,n 
c 

if (n . le . 0 ) return 

if (incx.eq.l.and.incy.eq.l)go to 20 
c 

c code for unequal increments or equal increments not equal 

c to 1 

c 

ix - 1 
iy - 1 

if (incx.lt. 0)ix - (-n+l)*incx + 1 
if (incy.lt. 0)iy = (-n+l)*incy + 1 
do 10 i ■ l,n 

stemp = c*sx(ix) + s*sy{iy) 
sy(iy) = c*sy (iy) - s*sx(ix) 
sx(ix) = stemp 
ix = ix + incx 
iy = iy + incy 
10 continue 
return 
c 

c code for both increments equal to 1 

c 

20 do 30 i * l,n 

stemp = c*sx(i) + s*sy(i) 
sy(i) * c*sy (i) - s*sx(i) 
sx(i) = stemp 
30 continue 
return 
end 

*** End of Subroutine SRot 


********************* *************************************************** 
subroutine srotg (sa f sb, c, s) 

******************************************* ***************************** 


c 

c 

c 

c 

c 


construct givens plane rotation, 
jack dongarra, linpack, 3/11/78. 

real sa, sb, c, s, roe, scale, r, z 

roe = sb 

if ( abs (sa) .gt. abs(sb) ) roe = sa 
scale = abs(sa) + abs (sb) 
if ( scale .ne. 0.0 ) go to 10 
c * 1.0 
s = 0.0 
r = 0 . 0 
go to 20 

10 r » scale*sqrt ( (sa/scale) **2 + (sb/scale) **2) 
r * sign (1.0, roe) *r 
c * sa/r 
s = sb/r 
20 z = 1.0 

if ( abs(sa) .gt. abs(sb) ) z = s 

if ( abs(sb) .ge. abs(sa) .and. c .ne. 0.0 ) z = 

sa = r 

sb = z 

return 


1.0/c 
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end 

*** End of Subroutine SRotG 


************************************************************************ 
subroutine sscal (n, sa, sx, incx) 

************************************************************************ 

c 

c scales a vector by a constant, 
c uses unrolled loops for increment equal to 1 . 

c jack dongarra, Unpack, 3/11/78. 

c 

real sa,sx(l) 

integer i, incx,m,mpl,n,nincx 
c 

if (n . le . 0 ) return 
if (incx.eq.l)go to 20 
c 

c code for increment not equal to 1 

c 

nincx = n*incx 
do 10 i * 1, nincx, incx 
sx(i) * sa*sx(i) 

10 continue 
return 
c 

c code for increment equal to 1 

c 

c 

c clean-up loop 

c 

20 m * mod (n, 5) 

if { m .eq. 0 ) go to 40 
do 30 i = l,m 
sx(i) = sa*sx(i) 

30 continue 

if ( n .It. 5 ) return 
40 mpl = m + 1 

do 50 i * mpl,n,5 


sx(i) = 

: sa*sx(i) 



sx (i + 

1) = sa*sx(i 

+ 

i) 

sx(i + 

2) = sa*sx(i 

+ 

2) 

sx (i + 

3) * sa*sx(i 

+ 

3) 

sx(i + 

4) = sa*sx(i 

+ 

4) 


50 continue 
return 
end 

*** End of Subroutine SSCal 


***************** ****************************************** ************* 
SUBROUT INE ZERO ( LEN , ARRAY ) 

************* ********************************************** ************* 

Just a little routine to zero the array. 

DIMENSION ARRAY (*) 

C 

DO 10 I- 1, LEN 
ARRAY (I) = 0. 

10 CONTINUE 
RETURN 
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END 

*** End of Subroutine Zero 


************************************************************************ 
SUBROUTINE TRIM (STRING, LSTRIN) 

************************************************************************ 

Return the length of STRING after trailing blanks, nulls, and tabs 
have been removed. 

CHARACTER* ( * ) STRING 

CHARACTER NULL, TAB 

Initialize the null and tab characters. 

NULL= CHAR < 0 ) 

TAB = CHAR (9) 

Loop backwards through the character string and find the last 
nonblank, nonnull character. 

LSTRIN= LEN (STRING) 

DO 10 L— LSTRIN, 1,-1 

IF (STRING (L : L) .NE.' ' .AND. STRING (L : L) .NE .NULL 
C .AND. STRING (L:L) .NE. TAB) THEN 

LSTRIN= L 
GOTO 20 
END IF 

10 CONTINUE 

ALL blank or null or tabs ! 

LSTRIN= 0 

20 CONTINUE 
RETURN 
END 

** End of Subroutine Trim 


*********************************************************************** 
SUBROUTINE UPCASE (STRING) 

******************************************************** *************** 
Convert this character string to upper case. 

CHARACTER* ( * ) STRING 
CHARACTER* 2 6 LOWER, UPPER 

DATA LOWER / 9 abcdef ghi jklmnopqrstuvwxyz' / 

C UPPER/ ' ABCDEFGHI JKLMNOPQRSTUVWXYZ ' / 

Don't worry about the trailing blanks - 

CALL TRIM (STRING, LSTRIN) 

Look for lower case letters and substitute upper case ones. 

DO 10 1= 1, LSTRIN 

LETTER= INDEX (LOWER, STRING (1:1) ) 


IF (LETTER. NE.O) STRING(I:I)= UPPER (LETTER: LETTER) 
10 CONTINUE 
RETURN 
END 

End of Subroutine UpCase 


*** 




